Why are my mesh points changing in my 1D linear convection C++ program?

Join the discussion
Ask a follow-up here, or get your own question answered by working scientists, mathematicians and engineers — people, not an autocomplete.
Real named experts · corrections over time · the nuance an AI answer skips
5 replies · 5K views
ssatonreb
Messages
5
Reaction score
0
I'm learning c++ and currently trying to write a little programme for 1D linear convection (wave equation).
I have managed set all boundary and initial conditions as well as a mesh.
But I have reached the point where I can't understand why the program does what it does.

In red I highlighted place where I am confused.

If someone could help me please.
In particularly I do not understand why does no matter what I do in the highlighted area of the code my mesh points are changed afterwards.

Code:
//  1D linear convection
#include <iostream>
#include <cmath>
using namespace std;

int main(){
    double nx=10;        // mesh points
    double nt=2;        // time steps 
    double dt=0.01;      // time step size
    double c=1;          // constant (wave propagation velocity)
    double dx=2/(nx-1);  // mesh size    
    int i,n ;
    
    const int NX = nx, NT = nt;     // array size assignment
    double x[NX], u[NT][NX];        // array declaration
    
    memset( u, 0, sizeof u );       // Zeroing out array u[NT][NX]
    memset( x, 0, sizeof x );       // Zeroing out array x[NX]
        
         // Construction of the mesh 
    for (i = 1 ; i < NX ; i++){ 
        x[i]=x[i-1]+dx;
         
    }
        // Initial and Boundary Conditions
    for (i = 0 ; i < NX ; i++){
        if ((x[i]>=0.5) && (x[i]<=1)){
                  u[0][i]=2;
        }
        else
            u[0][i]=1;
    }

// The problem lies here 
// --------------------------------------------------------------------------

    for (n=0; n<NT ; n++){
        for (i=1; i<NX-1; i++) {
            u[n+1][i]=u[n][i]-c*dt/dx*(u[n][i]-u[n][i-1]);
        }
    }
// --------------------------------------------------------------------------       
return 0;}

I really appreciate someones help.
 
Physics news on Phys.org
ssatonreb said:
I'm learning c++ and currently trying to write a little programme for 1D linear convection (wave equation).
I have managed set all boundary and initial conditions as well as a mesh.
But I have reached the point where I can't understand why the program does what it does.

In red I highlighted place where I am confused.

If someone could help me please.
In particularly I do not understand why does no matter what I do in the highlighted area of the code my mesh points are changed afterwards.
Change the line that I marked in red.
ssatonreb said:
Code:
//  1D linear convection
#include <iostream>
#include <cmath>
using namespace std;

int main(){
    double nx=10;        // mesh points
    double nt=2;        // time steps 
    double dt=0.01;      // time step size
    double c=1;          // constant (wave propagation velocity)
    double dx=2/(nx-1);  // mesh size[/color]    
    int i,n ;
    
    const int NX = nx, NT = nt;     // array size assignment
    double x[NX], u[NT][NX];        // array declaration
    
    memset( u, 0, sizeof u );       // Zeroing out array u[NT][NX]
    memset( x, 0, sizeof x );       // Zeroing out array x[NX]
        
         // Construction of the mesh 
    for (i = 1 ; i < NX ; i++){ 
        x[i]=x[i-1]+dx;
         
    }
        // Initial and Boundary Conditions
    for (i = 0 ; i < NX ; i++){
        if ((x[i]>=0.5) && (x[i]<=1)){
                  u[0][i]=2;
        }
        else
            u[0][i]=1;
    }

// The problem lies here 
// --------------------------------------------------------------------------

    for (n=0; n<NT ; n++){
        for (i=1; i<NX-1; i++) {
            u[n+1][i]=u[n][i]-c*dt/dx*(u[n][i]-u[n][i-1]);
        }
    }
// --------------------------------------------------------------------------       
return 0;}

I really appreciate someones help.

The line that I marked in red should be changed to this:
double dx=2.0/(nx-1); // mesh size

What was happening was that dx was being initialized to 0.0, which is not what you wanted. As you had it before, the expression on the right was 2/9, which is 0, due to integer division. That value was being promoted to 0.0 and stored in your dx variable.
 
Mark44 said:
Change the line that I marked in red.


The line that I marked in red should be changed to this:
double dx=2.0/(nx-1); // mesh size

What was happening was that dx was being initialized to 0.0, which is not what you wanted. As you had it before, the expression on the right was 2/9, which is 0, due to integer division. That value was being promoted to 0.0 and stored in your dx variable.

I think this is not a problem:
nx is a double, making (nx-1) a double, making 2/(nx-1) a double.

Another problem however, is that u[n+1] is indexing u out-of-range, which will give unpredictable results.
 
Thank you, Mark44.
I did as you suggested, but it didn't change my original problem, which is related with values of x.
Before I execute lines marked in red values of x are:
0 0.22(2) 0.44(4) 0.66(6) 0.88(8) 1.11(1) 1.33(3) 1.55(5) 1.77(7) 2
and after:
0 0.955 1 1.91203 1.99798 1.08797 1.00202 1 1 2

I have difficulty to understand the reason why does these values change at all because I don't do anything with these values.

Is it related with memory allocation?
 
Sorry - I glanced at nx = 10 and its comment and assumed it was an int. Since you are using it for the number of mesh points, why isn't it an int? When you use a variable to keep track of the number of things, use an integral type, like int or long.

I'm inclined to agree with I like Serena. It seems likely that your statement in red might be tromping over memory that isn't in the array. Your array u runs from u[0][0] through u[1][9].
In your nested loop the largest value of n is NT - 1 == 1 and the largest value of i is NX - 2 == 8.

The statement u[n+1] = ... attempts to store a value in u[2], which is outside this array.
 
I use double for nx because later I use it to calculate dx, which in case of int is just zero.
Thank you very much for your help!
I corrected my mistake and took into account Mark44's first comment I was able also get rid of a line:
const int NX = nx, NT = nt;
Also thank you, I like Serena.
I overlooked this point.
Thank you.