Anyone familiar with FDTD method?
- Thread starter TecHNooB
- Start date
Born2bwire
Diamond Member
- Oct 28, 2005
- 9,840
- 6
- 71
It is pretty straight forward. You take a differential equation and estimate using the discrete form of a derivative. When you have a multivariable differential equation, you can use a leapfrog method that will iterate the derivatives of one or more variables at a time (ie: space then time then space then time then...). There are considerations to be taken into account when in terms of the steppings for each variable to ensure stability of the algorithm.
For E&M, Yee's original paper is still pretty accessible and Taflove has a great book devoted entirely to FDTD. I even have my own 2D program I wrote in my spare time, the only thing I actually have on my Youtube account are some videos of the results that someone requested. You can really just write a FDTD code in an afternoon.
For E&M, Yee's original paper is still pretty accessible and Taflove has a great book devoted entirely to FDTD. I even have my own 2D program I wrote in my spare time, the only thing I actually have on my Youtube account are some videos of the results that someone requested. You can really just write a FDTD code in an afternoon.
Born2bwire
Diamond Member
- Oct 28, 2005
- 9,840
- 6
- 71
This is a response to a PM that I am putting in this thread because it may be useful to other people.
"Numerical solution of inital boundary value problems involving maxwell's equations in isotropic media" by Kane Yee: http://ieeexplore.ieee.org/sea...ta+%29&pos=1&access=no
I don't know if the above link will work but it is from IEEE Trans. Antennas and Propagation Vol 14, Issue 3, page 302 from 1966.
And a great reference is Taflove's book (Taflove coined the phrase FDTD I believe):
"Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd Edition.
http://www.artechhouse.com/Detail.aspx?strBookId=1123
Hell, just about any Computational EM book will have FDTD. Jin's FEM book, Weng Cho Chew's Inhomogeneous Media book (both are faculty for me so obviously I am biased but Taflove is THE book I think) have a brief sections on FDTD. Jin did FDTD early in his career so his book's section is not too bad but Chew really mentions it only for completeness I think.
All you do is take Maxwell's equations and express them in cartesian coordinates by each vector component. That gives you 6 differential equations with 6 unknowns. The currents are known, they are your excitation sources but you can also set any initial E and H fields that you want as an excitation. If you have a full 3D program then currents can replicate any desired excitation field but not so in a 2D program since the currents must be invariant in one dimension. Then you define your continuous parameter space (x,y,z,t) as being made up of discrete steps, a grid. You then use the basic approximations of first and second order derivatives to estimate the progression of the fields. The first trick is to understand the position of the fields in space on this grid. This is called the Yee grid and understanding this is key to FDTD. The E and H fields are not located at the same points in space. This is because it is more accurate to use a central difference approach to derivatives. So if we take two adjacent points to calculate the derivative, the resulting first-order derivative is located at the midpoint between these two points in central difference. Since the E field is dependent upon the derivatives of the H field, then the E fields will be located spacially between the H fields. Likewise in time, the E field will be calculated at points in time in between the times we calculate the H field. This dissassociation in time and space requires you to be careful when you define the bounday conditions and when you track the progression of the fields. However, it allows us to use leapfrogging to calculate the E and H fields separately.
The boundary conditions are the second problem. The easy thing to do is to have a PEC wall to encase your problem space. The problem here is that you now have made a cavity that will resonate (chances are it will resonate, of course there are situations where you it may not but most problem spaces will be many wavelengths in size and thus your frequency of choice will have a high chance of finding a mode.). Even if it isn't resonant, all the power will be reflected back into the problem space corrupting your solution. The easy way to solve this is to make the problem space large enough that by the time the first reflection comes back, you are already done with your measurement. This increases the problem size, increasing the memory and CPU costs and is undesirable.
The alternative is to use absorbing boundary conditions. You can use the Mur ABC conditions and/or Perfectly Matched Layer (PML). I prefer the PML because it does a better job than the Mur and it is easy to implement (a little tougher to implement it efficiently though). It requires you to split each field component into two extra components. You can do this over the whole problem space (easy but doubles the memory and CPU costs unnecessarily) or you can only do this at the PML (more work to code since you have to make a special set of arrays for the PML and figure out how to split and recombine the fields). Weng Cho Chew published a paper on how to implement the PML using coordinate stretching, this is easier than the original way of doing it. It just involves calculating some factors ahead of time to get the PML profile (just a few arrays) and then implementing these factors in the FDTD equations as coefficients. You do not have to worry about calculating extra derivatives at the boundaries and keeping old data like you do for Mur.
Stability is also an issue. You have to make sure that your time and spatial steps are small enough so that your solution will not blow up over time. Fortunately, the Yee algorithm already has the stability conditions defined for you but there will always be extra noise introduced by numerical error from floating point calculations and grid dispersion. These errors can make numerical noise that can eventually be blown up after enough iterations. So you may find you need to decrease the step sizes (try a factor of two or four) if you run into instability. Of course this means you have to increase the problems space and/or increase the number of iterations run to get the same simulation.
One tip about coding this, do not use multirank arrays for most programming languages. For example, in C++, a multirank array basically consists of arrays of pointers. A three rank array will contain a pointer to a pointer of an array. This requires many extra indirect dereferences and will quickly increase the amount of memory and CPU time required. Let's say I have a 512x512x512 problem space using a 64 bit machine and PML field splitting over the entire space. The arrays will be something like E[x][y][z][component]. You will need maybe 2 arrays like this. So that is what... umm... 16 GB just to hold all the pointers ([x] has 512 pointers pointing to [y] which has 512 pointers pointing to [z] which has 512 pointers to the arrays holding the vector components). Yeah, so inefficient multirank arrays will kill your speed and CPU. Do what I do, write an array class that will create and convert a multirank array interface to a single rank array. Then it is all hidden from you and you can tweak the array class as you like. In addition, split up the components into their own arrays too since a single rank array will occupy a contiguous block of memory (which probably isn't a problem for most problems that a class project code will do, you will usually hit the CPU time boundary on a PC than the memory boundary).
Writing the data output can take a lot of time (relatively) too. But in the end, you can use FDTD to do complex inhomogeneous medium problems (if you allow epsilon and mu to be part of the differential equations too) and can get these pretty time domain movie like this:
http://www.youtube.com/watch?v=GT1aRoauKm8&fmt=18
"Numerical solution of inital boundary value problems involving maxwell's equations in isotropic media" by Kane Yee: http://ieeexplore.ieee.org/sea...ta+%29&pos=1&access=no
I don't know if the above link will work but it is from IEEE Trans. Antennas and Propagation Vol 14, Issue 3, page 302 from 1966.
And a great reference is Taflove's book (Taflove coined the phrase FDTD I believe):
"Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd Edition.
http://www.artechhouse.com/Detail.aspx?strBookId=1123
Hell, just about any Computational EM book will have FDTD. Jin's FEM book, Weng Cho Chew's Inhomogeneous Media book (both are faculty for me so obviously I am biased but Taflove is THE book I think) have a brief sections on FDTD. Jin did FDTD early in his career so his book's section is not too bad but Chew really mentions it only for completeness I think.
All you do is take Maxwell's equations and express them in cartesian coordinates by each vector component. That gives you 6 differential equations with 6 unknowns. The currents are known, they are your excitation sources but you can also set any initial E and H fields that you want as an excitation. If you have a full 3D program then currents can replicate any desired excitation field but not so in a 2D program since the currents must be invariant in one dimension. Then you define your continuous parameter space (x,y,z,t) as being made up of discrete steps, a grid. You then use the basic approximations of first and second order derivatives to estimate the progression of the fields. The first trick is to understand the position of the fields in space on this grid. This is called the Yee grid and understanding this is key to FDTD. The E and H fields are not located at the same points in space. This is because it is more accurate to use a central difference approach to derivatives. So if we take two adjacent points to calculate the derivative, the resulting first-order derivative is located at the midpoint between these two points in central difference. Since the E field is dependent upon the derivatives of the H field, then the E fields will be located spacially between the H fields. Likewise in time, the E field will be calculated at points in time in between the times we calculate the H field. This dissassociation in time and space requires you to be careful when you define the bounday conditions and when you track the progression of the fields. However, it allows us to use leapfrogging to calculate the E and H fields separately.
The boundary conditions are the second problem. The easy thing to do is to have a PEC wall to encase your problem space. The problem here is that you now have made a cavity that will resonate (chances are it will resonate, of course there are situations where you it may not but most problem spaces will be many wavelengths in size and thus your frequency of choice will have a high chance of finding a mode.). Even if it isn't resonant, all the power will be reflected back into the problem space corrupting your solution. The easy way to solve this is to make the problem space large enough that by the time the first reflection comes back, you are already done with your measurement. This increases the problem size, increasing the memory and CPU costs and is undesirable.
The alternative is to use absorbing boundary conditions. You can use the Mur ABC conditions and/or Perfectly Matched Layer (PML). I prefer the PML because it does a better job than the Mur and it is easy to implement (a little tougher to implement it efficiently though). It requires you to split each field component into two extra components. You can do this over the whole problem space (easy but doubles the memory and CPU costs unnecessarily) or you can only do this at the PML (more work to code since you have to make a special set of arrays for the PML and figure out how to split and recombine the fields). Weng Cho Chew published a paper on how to implement the PML using coordinate stretching, this is easier than the original way of doing it. It just involves calculating some factors ahead of time to get the PML profile (just a few arrays) and then implementing these factors in the FDTD equations as coefficients. You do not have to worry about calculating extra derivatives at the boundaries and keeping old data like you do for Mur.
Stability is also an issue. You have to make sure that your time and spatial steps are small enough so that your solution will not blow up over time. Fortunately, the Yee algorithm already has the stability conditions defined for you but there will always be extra noise introduced by numerical error from floating point calculations and grid dispersion. These errors can make numerical noise that can eventually be blown up after enough iterations. So you may find you need to decrease the step sizes (try a factor of two or four) if you run into instability. Of course this means you have to increase the problems space and/or increase the number of iterations run to get the same simulation.
One tip about coding this, do not use multirank arrays for most programming languages. For example, in C++, a multirank array basically consists of arrays of pointers. A three rank array will contain a pointer to a pointer of an array. This requires many extra indirect dereferences and will quickly increase the amount of memory and CPU time required. Let's say I have a 512x512x512 problem space using a 64 bit machine and PML field splitting over the entire space. The arrays will be something like E[x][y][z][component]. You will need maybe 2 arrays like this. So that is what... umm... 16 GB just to hold all the pointers ([x] has 512 pointers pointing to [y] which has 512 pointers pointing to [z] which has 512 pointers to the arrays holding the vector components). Yeah, so inefficient multirank arrays will kill your speed and CPU. Do what I do, write an array class that will create and convert a multirank array interface to a single rank array. Then it is all hidden from you and you can tweak the array class as you like. In addition, split up the components into their own arrays too since a single rank array will occupy a contiguous block of memory (which probably isn't a problem for most problems that a class project code will do, you will usually hit the CPU time boundary on a PC than the memory boundary).
Writing the data output can take a lot of time (relatively) too. But in the end, you can use FDTD to do complex inhomogeneous medium problems (if you allow epsilon and mu to be part of the differential equations too) and can get these pretty time domain movie like this:
http://www.youtube.com/watch?v=GT1aRoauKm8&fmt=18
Born2bwire
Diamond Member
- Oct 28, 2005
- 9,840
- 6
- 71
Meh, my advisor once told me that you could write a Yee code in an afternoon and he is right. My 2D code that I wrote is around 750 lines, but most of that is just constructors, functions to read in a file for all the data, and functions to output the data. The actual FDTD portion of the code is only about 150 lines. The rest is just extra busy work I tacked on.
Finite difference codes in general tend to be short & simple. That's the appeal of the method... quick & dirty, but the applicability can be somewhat limited as things might get a little dicey with shocks, complex geometries, unstructured grids, etc.
<---is a finite element guy, lol
I mean besides having to sometimes do a fancy dance for boundary conditions, the basic core of finite difference methods is really simple: 1) decide on a finite difference stencil(s) for deriv approximation; 2) substitute into continuous equations; 3) write up PDE over the whole domain as a system of eqns.
A few other notes for this thread:
algorithm: it's often helpful with FD methods to draw out what the method is doing. That is, to find the solution at a point, you use some number of its neighbours to calculate spatial derivatives. Then this is linked through time by some time stepper (e.g. forward euler, RK-N, leapfrog). So if you draw out a picture in time & space, you can visualize how information is calculated at any given point & where that info is coming from.
stability: 1) Learn von-neumann analysis if you haven't already; it makes stability analysis for FD methods a cinch. But be warned: stability is NOT everything... classic example: crank nicholson (which is quite similar to the method being discussed here) is unconditionally stable, but it will produce spurious oscillations if the CFL number is chosen wrong.
2) leapfrog leaves much to be desired in terms of stability regions. iirc this is ok for equations arising from e&m, but be careful if you apply the same solver to other wave equation type things.
memory: you can gain significant speed-ups by coding your loops so that memory is accessed continuously (in addition to storing data continuously in memory). This is not always possible (e.g. you need to access A and A^T), but you'd be surprised how much difference reordering a few loops can make.
That said, I would always start by coding the algorithm in a way that makes the most sense to me. Once it works, start optimising.
I'm not that familar with E&M, but if solutions are smooth, it might be fun to experiment with a spectral spatial discretization and some kind of higher order time step like RK4.
- Sep 10, 2005
- 7,458
- 1
- 76
1 afternoon you say! I spent an entire weekend trying to unloack the mysteries of maxwell's equations. So far, all I've learned is that maxwell's equations can take on many ugly forms, PDE is much harder than ODE, and my vector calculus skills are pitiful. The rift between grads and undergrads is massive :/
MATLAB is perfect for FD implementation. Doing it in C shouldn't be too bad though - just make sure you're systematic about it.Originally posted by: TecHNooB
Hm.. finally starting to understand what's going on after working on the 1D case. However, keeping track of indices makes me sad
Originally posted by: CycloWizard
MATLAB is perfect for FD implementation. Doing it in C shouldn't be too bad though - just make sure you're systematic about it.Originally posted by: TecHNooB
Hm.. finally starting to understand what's going on after working on the 1D case. However, keeping track of indices makes me sad
Not if speed is an issue, no. Even if you're smart about it and take advantage of the compiled code for sparse matrices (i.e. not filling the matrix one element or even one row at a time), going to C can easily gain you a factor of 4-10 in performance.
That said, doing this stuff in matlab is WAY easier on the programmer. WAY easier.
- Sep 10, 2005
- 7,458
- 1
- 76
Originally posted by: eLiu
Originally posted by: CycloWizard
MATLAB is perfect for FD implementation. Doing it in C shouldn't be too bad though - just make sure you're systematic about it.Originally posted by: TecHNooB
Hm.. finally starting to understand what's going on after working on the 1D case. However, keeping track of indices makes me sad
Not if speed is an issue, no. Even if you're smart about it and take advantage of the compiled code for sparse matrices (i.e. not filling the matrix one element or even one row at a time), going to C can easily gain you a factor of 4-10 in performance.
That said, doing this stuff in matlab is WAY easier on the programmer. WAY easier.
I r doing it in matlab
Replicating this code in C would suck.
Born2bwire
Diamond Member
- Oct 28, 2005
- 9,840
- 6
- 71
Originally posted by: eLiu
Originally posted by: CycloWizard
MATLAB is perfect for FD implementation. Doing it in C shouldn't be too bad though - just make sure you're systematic about it.Originally posted by: TecHNooB
Hm.. finally starting to understand what's going on after working on the 1D case. However, keeping track of indices makes me sad
Not if speed is an issue, no. Even if you're smart about it and take advantage of the compiled code for sparse matrices (i.e. not filling the matrix one element or even one row at a time), going to C can easily gain you a factor of 4-10 in performance.
That said, doing this stuff in matlab is WAY easier on the programmer. WAY easier.
I have to agree. FDTD is very computationally intensive and there are few optimizations that you can do to speed it up. Since the program is very simple in the first place, it is usually better to write it in C, C++ or Fortran for speed. Matlab was much nicer for quick and dirty FEM and MOM codes. Matlab could generate a 2D mesh for my FEM problem, has a sparse matrix class to use, and it also will do the matrix solving. You can do this in C/++ and Fortran using the BLAS and LAPACK libraries but it is so much easier in Matlab.
Originally posted by: Born2bwire
Originally posted by: eLiu
Originally posted by: CycloWizard
MATLAB is perfect for FD implementation. Doing it in C shouldn't be too bad though - just make sure you're systematic about it.Originally posted by: TecHNooB
Hm.. finally starting to understand what's going on after working on the 1D case. However, keeping track of indices makes me sad
Not if speed is an issue, no. Even if you're smart about it and take advantage of the compiled code for sparse matrices (i.e. not filling the matrix one element or even one row at a time), going to C can easily gain you a factor of 4-10 in performance.
That said, doing this stuff in matlab is WAY easier on the programmer. WAY easier.
I have to agree. FDTD is very computationally intensive and there are few optimizations that you can do to speed it up. Since the program is very simple in the first place, it is usually better to write it in C, C++ or Fortran for speed. Matlab was much nicer for quick and dirty FEM and MOM codes. Matlab could generate a 2D mesh for my FEM problem, has a sparse matrix class to use, and it also will do the matrix solving. You can do this in C/++ and Fortran using the BLAS and LAPACK libraries but it is so much easier in Matlab.
Well, matlab can produce some nice, short FD code too:
spectral method for 1D Poisson's equation with period BCs, rhs() and n given:
h = 2*pi/n;
x = h*(1:n)';
fhat = fft( rhs(x) );
freq = [0:n/2, -n/2+1:-1]';
uhat = fhat./(freq.^2);
uhat(1) = 0;
u = real(ifft(uhat));
Or 1D convection-diffusion, homogenous dirichlet, 2nd order FD with rhs(), n, and diffusion coef mu given:
h = 2/(n+1);
x = (-1+h:h:1-h)';
K1D = spdiags(ones(n,1)*[-1*mu-h/2, 2*mu, -1*mu+h/2],-1:1,n,n)/h^2;
u = K1D\rhs(x);
Heh, anyway I should add that you'd get a 4-10x speed up over very well-written (read: vectorized to death) matlab. Over un-vectorized matlab, the speedup can be factors of 100 or even 1000.
The other big disadvantage of matlab is problem size. It just can't handle large matrices; clearly you can get farther by taking advantage of sparse() (if your problem is in fact sparse), but there's a limit... and for big problems, that limit is too low. Oh and there's "no" parallel matlab (there is star-p, but it isn't perfect and it isn't full-fledged).
Though you can greatly extend matlab's abilities with the MEX interface. MEX allows matlab to call compiled C code. I took a class in DG finite elements and while my matlab code was by far the fastest matlab, one student MEX-ed some parts of his code and was about 10x faster than I.
OP: if you ever get into doing scientific computing for a living... well, replicating that code in C will be life! Enjoy!
- Sep 10, 2005
- 7,458
- 1
- 76
Originally posted by: eLiu
Originally posted by: Born2bwire
Originally posted by: eLiu
Originally posted by: CycloWizard
MATLAB is perfect for FD implementation. Doing it in C shouldn't be too bad though - just make sure you're systematic about it.Originally posted by: TecHNooB
Hm.. finally starting to understand what's going on after working on the 1D case. However, keeping track of indices makes me sad
Not if speed is an issue, no. Even if you're smart about it and take advantage of the compiled code for sparse matrices (i.e. not filling the matrix one element or even one row at a time), going to C can easily gain you a factor of 4-10 in performance.
That said, doing this stuff in matlab is WAY easier on the programmer. WAY easier.
I have to agree. FDTD is very computationally intensive and there are few optimizations that you can do to speed it up. Since the program is very simple in the first place, it is usually better to write it in C, C++ or Fortran for speed. Matlab was much nicer for quick and dirty FEM and MOM codes. Matlab could generate a 2D mesh for my FEM problem, has a sparse matrix class to use, and it also will do the matrix solving. You can do this in C/++ and Fortran using the BLAS and LAPACK libraries but it is so much easier in Matlab.
Well, matlab can produce some nice, short FD code too:
spectral method for 1D Poisson's equation with period BCs, rhs() and n given:
h = 2*pi/n;
x = h*(1:n)';
fhat = fft( rhs(x) );
freq = [0:n/2, -n/2+1:-1]';
uhat = fhat./(freq.^2);
uhat(1) = 0;
u = real(ifft(uhat));
Or 1D convection-diffusion, homogenous dirichlet, 2nd order FD with rhs(), n, and diffusion coef mu given:
h = 2/(n+1);
x = (-1+h:h:1-h)';
K1D = spdiags(ones(n,1)*[-1*mu-h/2, 2*mu, -1*mu+h/2],-1:1,n,n)/h^2;
u = K1D\rhs(x);
Heh, anyway I should add that you'd get a 4-10x speed up over very well-written (read: vectorized to death) matlab. Over un-vectorized matlab, the speedup can be factors of 100 or even 1000.
The other big disadvantage of matlab is problem size. It just can't handle large matrices; clearly you can get farther by taking advantage of sparse() (if your problem is in fact sparse), but there's a limit... and for big problems, that limit is too low. Oh and there's "no" parallel matlab (there is star-p, but it isn't perfect and it isn't full-fledged).
Though you can greatly extend matlab's abilities with the MEX interface. MEX allows matlab to call compiled C code. I took a class in DG finite elements and while my matlab code was by far the fastest matlab, one student MEX-ed some parts of his code and was about 10x faster than I.
OP: if you ever get into doing scientific computing for a living... well, replicating that code in C will be life! Enjoy!
What the fork is that
And no, I'm going to build robots for the military
I've never seen a benchmark where ML is more than 2-3x slower than C code for any application, so I would really appreciate if you have any that indicate a greater slowdown. Now that I have a little more C under my belt, it might be worth my while to convert some of my more computationally-intensive programs if there's really that big of a difference for the type of work I'm doing.
Hm I guess I'm going off topic, but:
I've never seen any benchmarks, heh. All of my statements are from personal experience on problems arising from solving PDEs (e.g. finite difference, finite volume, finite element, and the more "how to solve Ax=b efficiently" things that go with it like CG, GMRES, multigrid, node reordering, preconditioners, etc).
Also, the conditions under which the C code is run matter *a lot*. Part of the advantage of using compiled code is all the magic the compiler can do for you. So if you compile -O0, I honestly have no idea how things will measure up. I've only ever tried -O2 or -O3 (not much difference between the two for my problems).
A few examples of matlab vs C:
1) I have a 1D finite volume code for solving the payne-witham traffic model (think shallow water eqns). It's as vectorized as possible (the only loop I have is for runge-kutta time stepping, which is unavoidable). C is 4-5x faster.
2) I think I mentioned this before, but I have a discontinuous galerkin code (2D problems like euler, navier stokes) that's pretty solid (I think). Someone else's partially MEX-ed code is about 10x faster.
3) Matrix assembly: for complex problems, this can be difficult to vectorize. And if you're having to visit every element of a matrix individually in matlab, the speed is going to be awful.
4) multigrid is easily 5x-10x slower (this is simply b/c the work in a basic multigrid routine is dominated by by smoothing... which usually entails using jacobi, gauss-seidel or SOR; these have no efficient implementation in m-code.
5) sparse matrix solves that take advantage of the matrix's sparsity pattern. For example, I know beforehand that my matrix is composed of 10x10 blocks; I can then re-write LU factorization to take advantage of this. The speed difference should be like 100x-1000x but this is a pretty unfair comparison.
Note that this is a -very- narrow class of problems. MATLAB and C have many more uses than solving PDEs (duh, lol) but that's where my experience lies.
I mean the fundamental problem is that matlab is interpreted. You lose twice in that respect: 1) no compiler optimizations for your specific code; 2) CPU is much less able to predict what will happen next. (Consider the simple example of a for loop to add the numbers in a big array. In C, the compiler will probably unroll that loop and maybe even call simd instructions. The CPU will be able to figure out that this loop is doing some kind of sequential access and prefetch data (among other things). In matlab, there's no compiler* and each loop iteration is separated from the next by the interpreter doing its thing.
Anyway I can believe that there are problems where C isn't more than 2-3x faster (some ideas discussed below).
Also keep in mind that beating MATLAB in C isn't necessarily trivial. I mean, suppose that your algorithm is in the form of a few nested for loops (think like... naive matrix-matrix multiply) and you code it like that in m-code. The same thing in C will be like 100x faster. Now instead consider the case where you've vectorized that m-code so that it's basically just A*B. Now matlab is calling the BLAS to do very efficient mat-mat mult and your nested for loops in C will have little or no advantage.
Or if your C code does silly things with memory allocation/access, uses data structures inappropriately, recurse like crazy, etc you could find that you won't gain much over matlab. Or if your algorithm is has certain kinds of bottlenecks, C might not win you too much (example: the matlab code consists of some simple problem setup and then calls to matlab's various compiled functions or your code has no significant loops). But even with all of that I'd be suprised if your code doesn't get at least a factor of 2 faster.
Last comments:
Coding for MATLAB is actually a lot like coding for a GPU. Speaking broadly, on a GPU, processing time is cheap and memory accesses are extremely expensive. Thus you want to do as much computation as possible for each visit to RAM. Similar idea in matlab: compiled code is very cheap and everything else is expensive. You want to do as many things as possible with the available compiled functions.
*Note: some of the examples I discussed above will be hard to reproduce because of the JIT (Just-In-Time) compiler that's new within the last version or two. In a handful of very specific situations, it can compile code on the fly to make things faster. So if you wanted to make a really simple for-loop speed comparison:
matlab:
sum = 0;
for i=1:100000
sum = sum + i;
end
C:
sum = 0;
for(i=1; i<=100000; i++)
sum += i;
The difference would be pretty slim b/c matlab will compile that loop. Instead you'd have to do something goofy like throw in a function call to evaluate (sum + i).
edit: don't get me wrong here. I love matlab and I use it all the time. Point is that the features that make matlab so convenient & easy to use are also the features that make it impractical for large problems.
Well if you remember my problem, it is basically a FD method. If Matlab/Maple really is a 5 to 10X more time I shudder to think the effect this would have on my code. Albeit poorly optimized it takes roughly 5 hours to do a run on a given set of parameters. Its also vaguely E&Mish, complex numbers, fun integrals etc.
Thanks for the reply. I'm not saying ML is good for all of these large problems - I've encountered my share of problems where it's very slow. I'm just trying to get a feel for the cases where it might be worth my while to do the programming in C rather than ML, given that it will take me approximately 10x as long (if not longer) to code in C.
What if I compile my ML? I am still running v7 and have the compiler...
lol, Maple will be like a billion times slower. For all of its awesomness at doing math with letters, it's -awful- at computationally intensive tasks. But it does have "arbitrary precision," which while not useful to me, is still kinda neat.
CycloWizard: Hmm well it's hard to make generalizing statements about when C will beat matlab significantly. What kind of problems do you work on? For virtually everything associated with solving PDEs, C will win. So in my case, unless I'm doing a one-off code or I'm just solving some small problems to demonstrate that an idea works, I'll ultimately go to C. (Since I often fail at C, my first step usually involves doing something in matlab so I can compare results. lol)
I mean I'd guess that the biggest underlying causes are loops and array accesses (which is just... more loops). Unless each loop is doing an epic amount of work per iteration in a minimal number of compiled function calls, I'd expect C to be a lot faster. By epic, I mean you're doing 10 iterations and each iteration takes 60 sec, as opposed to 100000 iterations and each iteration takes 0.006sec. And by 'minimum number of blah' I mean that each loop would consist of like... a few vector-vector or matrix-vector multiplies, a call to backslash and little more. (i.e. if I can do something with 10 matrix-vector multiplies or lump it together and just do 1 matrix-vector multiply, matlab is usually faster with the latter.)
That also ignores the time needed to set up the problem. This often happens one row/col or even one element at a time, and can be a substantial time-sink depending on the problem. But it's only significant if the "solving" part is simple (like no time stepping).
Also, what compiler are you talking about? You mean like EML? Or the (related) thing attached to simulink?
The problem that comes to mind is, ironically enough, a finite difference problem: four coupled PDEs for reaction-diffusion over long time scales. The goal was to determine the best-fit reaction-diffusion model parameters (I think there were nine?), so I had to solve the whole thing many, many times. Each time took a few hours I think, so I would just let it run for days on my home computer. The code was pretty simple and (I think) fully vectorized, but it was just a lot of computation. There was probably a more efficient method for solving the PDE's (like FE's
), which is probably where I could make the biggest difference, but when I was setting it up I didn't know much about other PDE methods, let alone coding anything by ML.The "MATLAB Compiler":
Doesn't seem to make a big difference in computation time, but that could be because I don't know anything about optimizing compilers so I just use default options for everything.
Originally posted by: CycloWizard
The problem that comes to mind is, ironically enough, a finite difference problem: four coupled PDEs for reaction-diffusion over long time scales. The goal was to determine the best-fit reaction-diffusion model parameters (I think there were nine?), so I had to solve the whole thing many, many times. Each time took a few hours I think, so I would just let it run for days on my home computer. The code was pretty simple and (I think) fully vectorized, but it was just a lot of computation. There was probably a more efficient method for solving the PDE's (like FE's
The "MATLAB Compiler":
Doesn't seem to make a big difference in computation time, but that could be because I don't know anything about optimizing compilers so I just use default options for everything.
Long time scales = how many iterations and what time stepping method? Can you estimate roughly how long each iteration takes and/or tell me what the size of the problem is (dimension of your matrix). Is the source term for reactions expensive to evaluate? Is your domain simple or do you employ mapping, embedded boundaries, or any other fancy dances?
Also, mcc is not quite what you think it is. Unless i'm wrong, mcc enables the running of m-code outside of matlab. It does this by making aspects of the matlab environment available outside of matlab through standalone libraries (i think this is called the MCR? I forget). I also remember some stuff happening with CTF files (some kind of encrypted thing) but I dunno wth was going on there. My understanding is that you probably shouldn't expect too big of a performance increase b/c the matlab core is still basically running in the background. Like when you start up an application created using mcc, I would expect that it'd take as long as starting matlab. Also I'd expect the memory usage to be comparable to running matlab (you can often do a lot better w/memory usage than matlab does.)
That said, I'm not super familiar with mcc, so I apologize if anything I've said is wrong. I tried using experimenting with it briefly to compare the performance of the code it generated with the C code that I wrote. Something didn't work so I moved on to more productive things, lol. Maybe I should revisit that. Anyway the only bit of advice I can offer is to run mcc with the -nojvm option. You won't be able to plot or do anything graphical, but that might speed the code up.
Pretty simple domain - just 1-d with one Neumann and one Dirichlet BC. 24*60*60*10 time steps (i.e. 0.1 s stepping for 24 hours) with only about 20 spatial nodes. The source term for this set of reactions was probably pretty expensive (A^a+b^b+C^c-D^d or something like that).
Yeah, that makes sense. I've only recently started compiling and writing GUIs for ML because I've started collaborating with some biologists. I'm used to just running from a script file, but I found out I could install MCR on their computers and give them something pretty to look at, which has made them pretty happy. From what I can tell, the CTF file is basically a zip file containing the auxiliary install files that are called by the executable. Since everything I've compiled is basically image processing, I give them a picture on the GUI to show what is going on (they're not too worried about performance I don't think).
20 spatial nodes? That's a tiny problem, lol. Is this implicit or explicit stepping? In either case, with nearly 1 million time iterations, I'd bet that you'll see a 4x speedup with C, maybe more. But if it's implicit and you solve a 20x20 system each time, you'll want to call LAPACK's tridiagonal or banded solver depending on what stencil you use. Take care that you only call pow() for fractional powers (except 1/2, then just use sqrt). For integer powers, compute by multiplication.
Oh ok. I thought you were trying to use mcc to make your matlab code go faster; my bad. haha... generating pretty pictures is half the point of doing computation
I guess I spend too much time with engineering types; I think everyone I know has a copy of matlab (maybe not legally but still).Yes, very tiny, but a lot of different time scales at work - fast, reversible reactions and slow to fast diffusion. I wrote explicit code because I didn't have much time to work on it (was originally a simple problem for a course project) and just let it run. And run. And run. I really had to shrink the time step a lot to get it to match up with the exact solution of the simplified case. If I ever get time to work on that particular problem again, I probably will go with C just for practice if nothing else.
I've found that one afternoon coding for a biologist will get me a paper, whereas solving a new engineering problem will take months. It's always fun to be able to solve someone else's huge problem with a few lines of code, saving them months of tracing things in ImageJ.Oh ok. I thought you were trying to use mcc to make your matlab code go faster; my bad. haha... generating pretty pictures is half the point of doing computation
- Sep 10, 2005
- 7,458
- 1
- 76
Hm. My 2-d fdtd code for maxwell's equations freaks out after several iterations (everything goes to infinity or -infinity). Any idea why this is happening? Just a prod in the right direction will do
I set all the boundaries to E = 0 and set one of the E vectors to a gaussian pulse. The matrices fill up and then blow up lol
Edit: actually, i think it might be a stability criterion. didnt put those in yet. ~_~
Edit2: still blows up
Edit3: Hm.. my magnitudes seem way to high, i wonder why that is O_O
I set all the boundaries to E = 0 and set one of the E vectors to a gaussian pulse. The matrices fill up and then blow up lol
Edit: actually, i think it might be a stability criterion. didnt put those in yet. ~_~
Edit2: still blows up
Edit3: Hm.. my magnitudes seem way to high, i wonder why that is O_O
TRENDING THREADS
-
Discussion Zen 5 Speculation (EPYC Turin and Strix Point/Granite Ridge - Ryzen 9000)- Started by DisEnchantment
- Replies: 24K
-
Discussion Intel Meteor, Arrow, Lunar & Panther Lakes Discussion Threads
- Started by Tigerick
- Replies: 19K
-
-
-
Facebook
Twitter