I have an iterative solver where it solves a PDE by given an initial guess and updating it until a convergence criteria is met. My question is embarrassingly simple. There's a point in the code where I update my old solution with the new for every iteration and this is done in MATLAB simply by: Uold = Unew; . However, this doesn't work in my C++ code and it breaks my code completely.

An example for the solver in C++:

This is an incomplete code, but the point I wanna show is that this line where I update my values is incorrect: //uoldk = unewh; //INCORRECT , basically I get nans . Where both of them are Eigen matrices defined as:
Thanks

basically I get nans

My first thought is to use a debugger to check for division by zero, or something close enough to zero to be a problem. It usually a good idea to explicitly code a check on the divisor to see it is not zero before doing the division. The same concepts apply to other functions that might bomb because the arguments are out of range or domain.

Edit:
Using the debugger, does it bomb on the very first iteration, or on one of the other values after that?

Do you validate your data? Edit: It may be quicker to write code to validate the data, rather than tediously debugging your way backwards through the code. Validation is really important, one should always do it.

Some other things that may help, having looked at the code from your other topic:

* Make use of version control programs like git. One can either use it from the shell, or there is integration with IDE's. I say this because your code had a lot of stuff commented out;

* If you have functions with lots of parameters, put each one on it's own line, rather than one giant line of code, as in.



This helps with display, especially on this site, the window for the code is fairly narrow. Also, if I am looking at the code in shell, it's easier if it doesn't word wrap.

* really try not to have giant functions (including main) it can be a real nightmare trying to break them up after the fact. Remember a function should only do one thing, if it's more than say 20 LOC, it's getting too big. I you are worried about passing lot's of variables about, consider putting them and their functions into a class or struct, so the functions can have direct access to the member variables
There are advantages in doing this. An object doesn't exist until it's constructor is completed. I like the class to inherit (with protected) a Validation class, whose functions I call from the constructor of the class in question. One can alos publicly inherit an Interface class, which has all the interface function in it and avoids cluttering the object's class.One can also throw exceptions. The downside of all this is that you may not want to have to learn all this other stuff;

* put a few lines of comments says what a function does;

* put a comment on the declaration of variables to say what they are. This is really helpful for others;

* if there is a link to documentation that explains what process you are following (wiki say) put it in a comment;

* Try not to use variables i and j or anything else that looks too similar. It will bite you one day when you get the order wrong. I prefer to use row and col if that is what they are;

* Do try to put your own code into it's own namespace. This will also save you when you inadvertently use an identifier that already exists in a library that you use, including the STL;

* one can use a namespace alias, as in this for example:

namespace ed = Eigen::Dense;

then qualify anything in Dense namespace with ed::
This allows one to use whatever abbreviation that suits for a namespace. This is a very common practice in the boost library;

I have an interest in this code, I am a Engineering and Mining Surveyor, and I am applying for a job doing data processing for geophysical surveys (Magnetic, Gravimetric and Spectral). It will be good for me to understand some of what goes on behind the scenes in the code. So I was wondering what application/purpose your code has in the spectral realm? I have been doing some study - mainly ruminating the wiki pages.

Regards :+)

Oh, and one more really big thing:

When compiling, always specify the c++ standard, and use -Wall -Wextra as a minimum. The compiler is your friend, use it to tell you about the problems.

It's worh reading the manual, there are a zillion options, but there are some useful ones that aren't included in Wall or Wextra.

but the point I wanna show is that this line where I update my values is incorrect: //uoldk = unewh; //INCORRECT , basically I get nans .

Why do you suppose the error is on that line?

What is in unewh before you try to assign it to something else?


BTW, where do you set phikmax?

You aren't really showing enough code or giving enough information about what you are trying to do to get meaningful answers.

I don't know what compiler/IDE is being used, it should be noted that Visual Studio 2022 integrates working with git/github into the IDE. Making it easier to do version control from the IDE.

I am still learning how to work with git and github so I really can't give pointers how to use the VS feature.

And just my opinion, using the 3rd party Eigen library isn't a beginner's topic. This is more general C++ programming.

@TheIdeasMan

My first thought is to use a debugger to check for division by zero, or something close enough to zero to be a problem. It usually a good idea to explicitly code a check on the divisor to see it is not zero before doing the division. The same concepts apply to other functions that might bomb because the arguments are out of range or domain. Edit: Using the debugger, does it bomb on the very first iteration, or on one of the other values after that?

Okay, I guess I should mention, that I have another code that I wrote in MATLAB that returns correct results and I am using it as a way to benchmark against if that makes sense. The reason, I suspect this line in particular is if I comment out this exact line from both codes I get exact results macthing from MATLAB and C++ code. Granted, the results when commenting this line are not correct sense the code will reach iterations=500 in MATLAB or count=501 in C++. This is why I believed the error is stemming from this particular line.

Do you validate your data? Edit: It may be quicker to write code to validate the data, rather than tediously debugging your way backwards through the code. Validation is really important, one should always do it.

I guess I do, as I mentioned I have a MATLAB code with a simple example that was tested and have been validated and verified as well and I use it to compare my C++ code with.

Also, you give excellent feedback on how to improve my code, believe me everything you recommended is my end goal here. But, I am VERY slow in learning everything about C++ and how to improve the efficiency. One of my big goals, is parallelization of the code and what you suggested will be VITAL to me. Thanks, I will slowly work my way to fix my aapx function.

Note: I do have a remote repo where I keep many copies of my code safe.

I have an interest in this code, I am a Engineering and Mining Surveyor, and I am applying for a job doing data processing for geophysical surveys (Magnetic, Gravimetric and Spectral). It will be good for me to understand some of what goes on behind the scenes in the code. So I was wondering what application/purpose your code has in the spectral realm? I have been doing some study - mainly ruminating the wiki pages.

Thanks! I should clarify, that "spectral" here is referred to as spectral methods which are numerical methods in solving PDEs. My code has space weather applications and is meant to study ionospheric turbulence and their impact on communication systems such as GPS.

@George P

And just my opinion, using the 3rd party Eigen library isn't a beginner's topic. This is more general C++ programming.

I was not sure, I felt pretty lost and that my question was basic, but you're probably right. Thanks.

Do not expect matlab and C++ to match exactly. There will be small (insignificant) variation in the results due to little things like order of operations or exact algorithm used and so on.

That said, checking it against it at the same step is the way to go. But I don't see the problem yet.

@jonnin

That said, checking it against it at the same step is the way to go. But I don't see the problem yet.

Actually that's what I did, and I found that the C++ and MATLAB solution are the same but with a factor of 1000 off. I guess, I will be printing everything out and finding where that factor is coming from. Thanks.

JamieAl wrote:
I guess I do, as I mentioned I have a MATLAB code with a simple example that was tested and have been validated and verified as well and I use it to compare my C++ code with.

JamieAl wrote:
I guess, I will be printing everything out and finding where that factor is coming from. Thanks.

I meant validate as in writing C++ code to check that the input is in the correct range.

Printing out values is called "Poor man's debugging", there are a bunch of reasons why this becomes painful: time taken to write code to print values; time taken to remove that code later are just two of them. Consider using an actual debugger to step through the code, looking at how values change, and deduce where it all goes wrong.

Using Version Control (VC) is a different thing to storing multiple copies elsewhere. With VC one checks out the latest version of the code, do work on it, then check it back in again with a comment about what you did. It's possible to revert to previous versions. In a team environment it's possible to branch and merge. Consider looking at git, it's one of the best out there. And it's easy to use inside an IDE as George P mentioned.

There are some other things I noticed:

counter = counter + 1;

can be written as

counter++;

With std::pow, it tends to be terribly inefficient for integer exponents, maybe (I am guessing, it's not set in stone how it is implemented) because it does a binomial series to evaluate it. If you are just squaring, then simply multiply the number by itself:



Note I had 2.0 instead of 2, it saves the compiler from doing an implicit cast, and reinforces the idea we are using doubles.

@TheIdeasMan
Thanks, I appreciate it! I will fix all of this once I am done with setting up my functions.

@JamieAl

Hi,

Just wondering what a PDE is?

If you post code for your current problem, we can help you better.

In the file Test.cpp

There are a number of places where you could make things constexpr:

constexpr double Lx = 2*EIGEN_PI; //2*M_PI;

constexpr double dx = Lx/nx; // nx is staic const

If dx/2.0 is used a lot make it constexpr too:

constexpr double halfdx = dx/2.0;

I find 2.0 easier to read than 2.

With this:

If you set your editor to convert tabs to spaces (tab is 4 or less spaces) it will display a lot better here. On this site a tab is 8 spaces.

Lx is 2 * PI, so

KX[i] =2.*EIGEN_PI*i/Lx; //2.*M_PI*i/Lx;
becomes

KX[i] =i; //2.*M_PI*i/Lx; Lx and 2 * M_PI cancel out

similar with
KX[i] = 2.* EIGEN_PI * (-i + 2.*k_counter) / Lx;
becomes
KX[i] = -i + 2.*k_counter ;

This:



Can be written on one line:

double* KX = (double*) fftw_malloc(nx*nyk*sizeof(double));

There are a lot of similar instances.

I hope this helps.

I discovered what a PDE is Partial Differential Equation. Sorry for being a goose, all I had to do was Google it!! :+)

@TheIdeasMan

I discovered what a PDE is Partial Differential Equation. Sorry for being a goose, all I had to do was Google it!! :+)

np! I should've been more clear and also thanks for all the help again. I probably will be back for more questions or to get more suggestions like this for the sake of efficiency of my code.