I'm writing a program to solve the Bloch equations for a magnetic field that changes in time. The values of my magnetic field vary over too great a range to use a simple Runge Kutta method, but I have a standard 4th order with a constant step size set up to work with a constant uniform magnetic field. Now I'm trying to rig an a program that uses an embedded 4-5 order method to use the truncation error, and adjust my stepsize (ala the methods in Numerical Recipes in C (2nd ed))
The function as I have it is thus:
void rungekuttaKC(double a[6], double b[6][6], double c[6], double cs[6], double h, double gammaN, double bx, double by, double bz, double *mu) {
double *mu5 = new double[3];
mu5 = mu;
double *mu4 = new double[4];
mu4 = mu;
double *kx = new double[6];
double *ky = new double[6];
double *kz = new double[6];
double add[3] = { 0.0, 0.0, 0.0 };
double addmu[3] = { 0.0, 0.0, 0.0 };
double addmus[3] = { 0.0, 0.0, 0.0 };
for (int i = 0; i < 6; i++) {
for (int j = 0; j < i; j++) {
add[0] = add[0] + kx[j] * b[i][j];
add[1] = add[1] + ky[j] * b[i][j];
add[2] = add[2] + kz[j] * b[i][j];
}
kx[i] = h*gammaN*((mu[1] + h*add[1]) * bz - (mu[2] + h*add[2]) * by);
ky[i] = h*gammaN*((mu[2] + h*add[2]) * bx - (mu[0] + h*add[0]) * bz);
kz[i] = h*gammaN*((mu[0] + h*add[0]) * by - (mu[1] + h*add[1]) * bx);
}
for (int l = 0;l < 6;l++){
addmu[0] = addmu[0] + kx[l] * c[l];
addmu[1] = addmu[1] + ky[l] * c[l];
addmu[2] = addmu[2] + kz[l] * c[l];
addmus[0] = addmu[0] + kx[l] * cs[l];
addmus[1] = addmu[1] + ky[l] * cs[l];
addmus[2] = addmu[2] + kz[l] * cs[l];
}
for (int m = 0;m <=2; m++) {
mu4[m] = mu4[m] + addmus[m];
mu5[m] = mu5[m] + addmu[m];
mu[m] = mu4[m];
}
//delete[] mu5;
//delete[] mu4;
//delete[] kx;
//delete[] ky;
//delete[] kz;
//delete[] add;
//delete[] addmu;
//delete[] addmus;
}
a, b, c, and cs are the values in the Bucher tableau for the Cash-Karp method, h is the step size,1.0e-8.
static double a[6] = { 1.0 / 5.0, 3.0 / 10.0, 3.0 / 5.0, 1.0, 7.0 / 8.0 };
static double bij[6][6] = { {0.0,0.0,0.0,0.0,0.0,0.0},
{1.0 / 5.0, 0.0, 0.0, 0.0, 0.0, 0.0},
{3.0 / 40.0, 9.0/40.0, 0.0, 0.0, 0.0, 0.0},
{3.0/10.0, -9.0/10.0, 6.0/5.0, 0.0, 0.0, 0.0},
{-11.0/54.0, 5.0/2.0, -70.0/27.0, 35.0/27.0, 0.0},
{1631.0/55296.0, 175.0/512.0, 575.0/13824.0, 44275.0/110592.0, 253.0/4096.0} };
static double c[6] = { 37.0 / 378.0, 0.0, 250.0/621.0, 125.0/594.0, 0.0, 512.0/1771.0 };
static double cs[6] = {2825.0/27648.0, 0.0, 18575.0/48384.0, 13525.0/55296.0, 277.0/14336.0, 1.0/4.0 };
the magnetic field values are set up outside the function by an interpolator, and passed in, but just to make sure it's producing oscillations, I am testing it with a constant field, where bx and by are 0, and bz is 5. The mu that gets passed in starts as {1,0,0}
Now, the whole program runs in the debugger when I don't use any of the delete commands at the end of the function, but there's a pretty obvious memory leak.
When I use delete with only the arrays I initialize with new, a window pops up with the error "Main.exe has triggered a breakpoint." and no other information. In the output window I get " HEAP[Main.exe]: Invalid address specified to RtlValidateHeap( 00470000, 0016F514 ) Main.exe has triggered a breakpoint."
And when I delete all of them at the end of the function, I get the same error, but a different address in the RtlValidateHeap.
I've tried running in release mode, and switching my runtime library and I get "Critical error detected c0000374 Main.exe has triggered a breakpoint."
Am I misusing delete[] or is there something else I'm doing wrong? what does this error mean?
You are deleting add, addmu, and addmus. These were allocated on the stack (not on the heap with new). Therefore, they cannot be delete[]'d.
Since you know all the sizes of your arrays at compile-time and they are small, you should probably allocate them all on the stack. For dynamic arrays, in C++ the best practice is to use a container like vector, which frees you from the responsibility of using delete manually.
Additionally, as Borgleader said, the lines
double *mu5 = new double[3];
mu5 = mu;
double *mu4 = new double[4];
mu4 = mu;
are not doing what you expect. First, a double[3] is allocated on the heap, and its address is given to mu5. Then, you erase that address and overwrite it with mu. When you delete[] mu5, you are freeing the memory pointed to by mu (which is probably not intended), and the new double[3] is still standing. To correct it now with a minimum of changes, you could replace that snippet with
double *mu5 = new double[3];
std::copy(mu, mu + 3, mu5);
double *mu4 = new double[4];
std::copy(mu, mu + 4, mu4);
However, writing this is C style, not C++ style (although it uses C++ syntax).
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With