How can I limit the array in MathCAD to prevent crashing?

[haitrungdo82]

haitrungdo82: I agree with all of your statements in all of your posts yesterday, except for the 200 min. An analysis time of 30 min is better. Your wind velocity is 4.47 m/s. If you run the simulation more than, say, 30 min, the pollutant has already been moved off of the 10 km analysis space by the wind. Therefore, it is irrelevant to run the simulation longer.

D = 0.000 017 (m^2)/min appears incorrect. A value this small would have no effect on the problem, and would be irrelevant.

First, just for fun, try the following.

(1) nt := 2000
(2) nx := 200
(3) deltat := 6 s
(4) deltax := 50 m
(5) v := 0 m/s
(6) D := 0.01 (km^2)/min
(7) f_0,nx := 0
(8) f_0,0.5*nx := 1

Then try the following.

(9) nt := 2000
(10) nx := 2000
(11) deltat := 0.90 s
(12) deltax := 10 m
(13) v := 4.470 m/s
(14) D := 0.0017 (km^2)/min
(15) f_0,nx := 0
(16) f_0,0.5*nx := 1

You could also try the following, but it does not appear to be better than lines 9 through 16. Therefore, lines 9 through 16 seem better.

(17) nt := 3600
(18) nx := 2200
(19) deltat := 0.50 s
(20) deltax := 9 m
(21) v := 4.470 m/s
(22) D := 0.0017 (km^2)/min
(23) f_0,nx := 0
(24) f_0,0.5*nx := 1

Hi nvn,

Thumb up! All the negative values disappeared and the model is physically perfect! It seems like the diffusion coefficient is the cause of all messy, unphysical things. I'm just curious! What were you thinking when you were trying to eliminate the negative values? I'm very curious about this! Thank you!

 

[nvn]

haitrungdo82: My thinking was very vague. We know 0.5*(f_t,x+1 - f_t,x-1)/deltax is the first-order finite difference approximation of the first derivative, if I recall correctly. And v*deltat is an x (horizontal) distance. Therefore, this horizontal distance times the above-mentioned first derivative (which would be the tangent of the slope angle) should be a vertical (y) distance, if I recall correctly.

We know (f_t,x-1 - 2*f_t,x + f_t,x+1)/deltax^2 is the first-order finite central difference approximation of the second derivative (curvature), if I recall correctly. Therefore, similarly, when we multiply this value times D*deltat, we get another y distance. I do not recall the exact details.

Therefore, I vaguely envisioned that if we supply a D value too large or too small, the equation was overshooting below y = 0, or overshooting in the horizontal direction (?), which, in either case, would give garbage. I envisioned that if D is too small, it has no effect on the dispersion. If D is too large, the equation overshoots. I envisioned if deltax is too small, we get almost a singularity, called overflow. We want to increase D, but I envisioned that we cannot increase it too much unless we decrease deltat; otherwise, the problem will overshoot, as explained above. I envisioned that you cannot decrease deltat below a certain value, or else the problem will crash or overshoot. And I envisioned that you cannot decrease deltax below a certain value, or else the problem will crash or overshoot. E.g., if we want to simulate out to 10 km, it seems deltat cannot be less than 0.5 s, and deltax cannot be less than ~8 m. I do not understand very well what was going on, nor how to control the diffusion equation. If you uncover anything else, feel free to let us know.

I had to increase D seemingly too much (?); but that is the only way I could get the problem to run. In your given problem, D is overpowered by the relatively high wind velocity.

In the "Graph from Question" file, we see a lot of diffusion, and almost no wind. We can see that the centerline of the probability density moves slightly in the negative x direction, which means there is a very, very slight wind in the negative x direction.

By the way, to reply to a post, did you know you can just press the "New Reply" button? You don't need to press the "Quote" button on every reply. Also, always leave a space between a numeric value and its following unit symbol. E.g., 10 km, not 10km. See the international standard for writing units[/color] (ISO 31-0[/color]).

If you uncover anything else, feel free to let us know.

 

[haitrungdo82]

haitrungdo82: My thinking was very vague. We know 0.5*(f_t,x+1 - f_t,x-1)/deltax is the first-order finite difference approximation of the first derivative, if I recall correctly. And v*deltat is an x (horizontal) distance. Therefore, this horizontal distance times the above-mentioned first derivative (which would be the tangent of the slope angle) should be a vertical (y) distance, if I recall correctly.

We know (f_t,x-1 - 2*f_t,x + f_t,x+1)/deltax^2 is the first-order finite central difference approximation of the second derivative (curvature), if I recall correctly. Therefore, similarly, when we multiply this value times D*deltat, we get another y distance. I do not recall the exact details.

Therefore, I vaguely envisioned that if we supply a D value too large or too small, the equation was overshooting below y = 0, or overshooting in the horizontal direction (?), which, in either case, would give garbage. I envisioned that if D is too small, it has no effect on the dispersion. If D is too large, the equation overshoots. I envisioned if deltax is too small, we get almost a singularity, called overflow. We want to increase D, but I envisioned that we cannot increase it too much unless we decrease deltat; otherwise, the problem will overshoot, as explained above. I envisioned that you cannot decrease deltat below a certain value, or else the problem will crash or overshoot. And I envisioned that you cannot decrease deltax below a certain value, or else the problem will crash or overshoot. E.g., if we want to simulate out to 10 km, it seems deltat cannot be less than 0.5 s, and deltax cannot be less than ~8 m. I do not understand very well what was going on, nor how to control the diffusion equation. If you uncover anything else, feel free to let us know.

I had to increase D seemingly too much (?); but that is the only way I could get the problem to run. In your given problem, D is overpowered by the relatively high wind velocity.

In the "Graph from Question" file, we see a lot of diffusion, and almost no wind. We can see that the centerline of the probability density moves slightly in the negative x direction, which means there is a very, very slight wind in the negative x direction.

By the way, to reply to a post, did you know you can just press the "New Reply" button? You don't need to press the "Quote" button on every reply. Also, always leave a space between a numeric value and its following unit symbol. E.g., 10 km, not 10km. See the international standard for writing units[/color] (ISO 31-0[/color]).

If you uncover anything else, feel free to let us know.

Hi nvn,

Thank you very much for your advice! I was afraid that if I did not press the "Quote" button, you would not be aware of my new messages. If this is not the case, I will press the "New Reply" next time.

This is probably my last question for you. I was going to take the natural log (ln) of each value in the array. However, since there are many "0"s, the program will crash. There are 2 ways I can solve this problem. I can either convert all "0" into "1" or skip the "0" values. I don't know how to skip some selected values in the array, so I chose the first method.

I used "if" to convert all "0" to "1" as you can see from the attachment. However, only the upper part of the array was actually converted. The rest still contained "0". I copied 2 different areas of the same array in the attachment to illustrate my point.

So, my question is whether I missed anything when I applied the condition. Or if you know how to skip some values in an array, it would be great and I would like to learn.

Thank you very much for your help

 

[haitrungdo82]

You'll have to decide if negative values have proper physical meaning for what is being modeled. If not, there is some sort of problem, either with the model itself, the way it is set up (parameters), or the interpretation of the results.

If this is a well-known model, then there are probably papers written which discuss its use, stability, and overall behavior.

Hi gnell,

I just posted a question regarding how to convert some values in an array to other values or how to skip some values in an array. Hope you or nvn can help me since both of you are very good at MathCAD. Thank you very much!

 

[nvn]

haitrungdo82: We will indeed be aware of your new messages if you press "New Reply."

Try changing your program to something like f5_t,x := {f5_t,x if f5_t,x > 10^-8; 1 otherwise}, and see if it then works. You can adjust the -8 exponent to any value you prefer, such as -6, -8, -10, etc. Also, ensure you are displaying enough decimal places, or scientific notation, in your matrix to enable you to see the actual values stored in your matrix.