blockCoupled solver
I'm writing a solver that uses the blockCoupledLduMatrix, only the number of coupled variables is arbitrary  reactions are read from a text file.
Is there such a thing as a volVectorNField, where N is left arbitrary? Or do I have to make an explicit instantiation for a specific number? 
Hi Dave, Unfortunately the design of the block solver requires you to instantiate a specific type. You can make this a little easier on yourself by creating a template class for your system of equations. This way you can reuse it for an arbitrary number of equations. e.g. Let's say you create a class for your solver
Code:
template<class cmpt, int length> And then in your top level solver you can do something along the lines of... Code:
switch nEqns Regards Ivor 
Just to follow up a little more. The volVectorNField you ask about has been defined for a number of types. In the library VectorN in OpenFOAM1.6ext you will find definitions for volVector2Field, volVector4Field, volVector6Field and volVector8Field.

Hello you two,
I'm working with the block solver, too and found a workaround for that. Not that it is very shiny but for me it works. During run time a cannot change the number of species but during compile time. In my solver folder I added an additional header, which will do the definitions of the vectorN / tensorN. Before compile I just need to tell the template the number of species, the rest is done automatically. Code:
const int N=3; In my solver I just use the "mySomethingN", so at the toplevel it is the same for all number of species. I plan to write a script which will work as the executable and which will change the "N", compile the solver and then run the case. I just did not do it up to now, because I didn't mind recompiling on my own. Best Kathrin 
A follow up...
Should we just collect this discussion on the extendproject webpage? We could start a group there with special interest on the block solver! Best Kathrin 
Thanks Ivor and Kathrin.
Kathrin, what includes are you putting above that code? I'm getting an annoying compile error where it complains that it doesn't know what DiagTensorN and SphericalTensorN are... but I have included the definitions (I even compiled with the E option to see what the preprocessed code looked like, and sure enough, DiagTensorN is well defined before I try using it). I'm all for the blockCoupled SIG. 
Hi David,
here are the includes Code:
#ifndef myNThings_H The file holds other useful stuff but it is too big to post it here. I just opened such a working group on the extend portal and I can give you access to the file there! I invited you and Ivor to join! Best Kathrin 
David,
Have a look at the neutronDiffusion solver that I supplied at the OpenFOAM workshop. It uses Kathrin's approach and should work out the box with OF1.6ext. Regards Ivor 
Thanks for the suggestions, Ivor.
I'm trying two approaches and getting different errors for each. My system of equations (currently) involves 37 variables (but it changes). Is this too big? My first approach is based on Kathrin's suggestion. The error I get here is: Code:
error: variable ‘Foam::BlockLduMatrix<Foam::VectorN<double, 37> > blockM’ has initializer but incomplete type". I went to the VectorN library and tried to create an instantiation for 37 variables. I get "error: template instantiation depth exceeds maximum": Code:
.../src/OpenFOAM/lnInclude/VectorSpaceM.H: In static member function ‘static void VectorSpaceOps<N, I>::eqOp(V1&, const V2&, EqOp) [with V1 = Foam::VectorSpace<Foam::TensorN<double, 37>, double, 1369>, V2 = Foam::VectorSpace<Foam::TensorN<double, 37>, double, 1369>, EqOp = Foam::eqOp<double>, int N = 1369, int I = 193]’: 
37 variables!!! That's way bigger than I've ever tried.
You did the right thing trying to instantiate VectorN<scalar, 37>. You will probably have to do the same with TensorN<scalar, 37> and possible DiagTensorN and SphericalTensorN. This will depend on your problem type. I'm not surprised you received the 'exceeds maximum depth' message though. You will need to update the compiler options to get it to accept such a large number. In your Make/options file add the option ftemplatedepthXXX to the EXE_INC option. Here XXX will need to be set sufficiently large for the compiler to do it's thing. At an initial guess XXX would have to be at least 37^2. And yes, the compiler will be VERY slow with such a large number. 
Just some more food for thought. You might want to think seriously about whether you need the blockcoupled solver for this problem. Although in all likelihood it will work, with 37 coupled variables you're computational overheads could go through the roof. If this is a fully blockcoupled system then you run the risk of multiplying both your memory and computational requirements by a factor of 1000. I'm guessing the coefficients themselves will be relatively sparse so you might be wasting a lot of resources by storing 37x37 tensor coefficients. If your system is not very stiff you might want to be looking at other approaches like explicit RungeKutta and predictorcorrector type methods instead.
That being said, I'm also interested to see how the blockcoupled solver handles this one. Please keep us informed. 
I managed to compile VectorN for 2,4,6,8 and 37 variables. The libso file is 18 Mb. I will try writing the solver based on the blockCoupledScalarTransportFoam template.
How much calculation overhead? Would it be slower than directsolving the 37x37 reaction matrix on a cellbycell basis across the whole mesh 10 to 1000 times per timestep? My problem requires a pointimplicit solver only. I've calculated that using the blockCoupled solver requires a minimum of 16 kB per cell (plus it might cache the coefficients or record .oldTime() values or something). Since the reactions have a large length scale, I can use a coarse grid (say 1000 to 10,000 elements), so I'm not too concerned about memory. I'm more concerned about calculation time. I'm currently solving the reactions on a cellbycell basis. Transport is a problem as this requires the neighbours, so what I do is use the old value, and solve the whole mesh until neighbours match. This requires many iterations and relaxation. I've tried using lduMatrix but it failed. The reaction matrix is sparse (not as sparse as a transport lduMatrix, but close), square, but not diagonally dominant. Even if I enforce diagonal dominance, it is too unstable. These reactions are stiff. I also tried using an ODEtype of solution, but the Courrant number requirement (Co < 1) limits the timestep too much, so it would take years to achieve a solution. 
It seems this might be your best bet then.
In terms of calculational time, the matrix preconditioners (both diagonal and Cholesky) require inversion of the matrix coefficients for each solution iteration. Since you're talking about pointimplicit coupling, this will be roughly two 37x37 tensor inversions for each cell for each matrix iteration (assuming your matrix is asymmetric). If your iterative solver is going to take 1000 iterations to converge then, yes, this could be as slow as the explicit approach. Based on what you've said, though, I think you could see major improvements over your previous efforts. 
I don't want to hijack this thread, but are there any examples of the blockcoupled solver used with multiple regions? All the folks using openFoam for electrochemical simulations (batteries and fuel cells) are dying for this (ex. we have 3 regions and ~3 coupled pde's spanning all 3 regions).

No specific examples yet, but the underlying framework is there. I've done something along these lines for the simplified P3 equations and might be able to give some guidance. Are the equations the same over all three regions, i.e. do you need to split the problem into three regions or can you just change the material properties. If you're really interested in this, post a new thread and we can start separate discussions.

nonsearchable discussion
Hi,
Why move away from an open, Googlesearchable discussion to a closed, nonsearchable discussion? Don't get me wrong, I like the extendproject, but I don't see the benefit of having an additinal forum/discussionsite which isn't even searchable. regards, Sylvester Quote:

Hi Sylvester,
You make a good point. I have mixed feelings on the topic. First, I agree the discussion we have on this should be accessible. These forums are very accessible, and the Extend Portal has some problems:  Interface issues (I find it difficult to navigate the portal)  Nonsearchable (alleged) I'm not sure why the portal isn't searchable. It should be possible to do, given that other database / phpinterface websites are (e.g. wikipedia). If the Portal wants to displace the cfdonline forums, users have to find it easy to switch to. These issues won't help. On the other hand, as you may be aware, there are two variations of OpenFOAM: Official and Extend. It is unfortunate that these two groups don't get along very well. I'm not sure how this will play out, but there's a chance that the groups will split as we've seen in the past with other similar situations (e.g. MySQL and Drizzle). If that happens, the Portal is the likely candidate of where Extendrelated discussions will be taking place, so I think it is important to participate. Participating in the portal also includes providing constructive feedback on its interface / function, in an effort to improve it for everyone. 
Quote:
Ah. And look for the last posting by "official" (apart from the announcement of 2.0). Seems like some people already gave up on the Forum (in favour of Twitter) 
Fair enough. The "Communitydriven Collaboration" section of the portal looks very functional too.

Back on topic, I just finished uploading a slightly more generalized implementation of the VectorN library. The files are on the Portal... somewhere... arrrggg... is there a hard link to the "Block Solver in OpenFOAM" project I can put on here?
Let's see... go to the "Communitydriven Collaborative Projects" link on the side bar under community. Then choose "Block Solver in OpenFOAM". Files are available from a link on the top. 
All times are GMT 4. The time now is 16:47. 