# blockCoupled solver

 Register Blogs Members List Search Today's Posts Mark Forums Read

 July 6, 2011, 17:34 blockCoupled solver #1 Senior Member   David Gaden Join Date: Apr 2009 Location: Winnipeg, Canada Posts: 436 Rep Power: 15 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?

 July 6, 2011, 21:44 #2 Member   Ivor Clifford Join Date: Mar 2009 Location: Switzerland Posts: 91 Rep Power: 10 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 mySystemOfEquations { typedef VectorN vectorType; ... blockSolverPerformance solve(); }``` You can use the length template parameter to do looping, etc. or whatever you need to construct the block matrix. And then in your top level solver you can do something along the lines of... Code: ```switch nEqns { case(2): mySystemOfEquations eqns(...); eqns.solve(); case(4): mySystemOfEquations eqns(...); eqns.solve();``` And so forth. I hope this helps you. There may be more elegant approaches, like using runtime selection but the above approach is what I've used in the past. Regards Ivor mm.abdollahzadeh likes this.

 July 6, 2011, 21:48 #3 Member   Ivor Clifford Join Date: Mar 2009 Location: Switzerland Posts: 91 Rep Power: 10 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 OpenFOAM-1.6-ext you will find definitions for volVector2Field, volVector4Field, volVector6Field and volVector8Field.

 July 7, 2011, 03:33 #4 Senior Member   Kathrin Kissling Join Date: Mar 2009 Location: Besigheim, Germany Posts: 134 Rep Power: 10 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; typedef VectorN myVectorN; typedef TensorN myTensorN; typedef DiagTensorN myDiagTensorN; typedef SphericalTensorN mySphericalTensorN; // typedef symmTensorN mySymmTensorN; typedef dimensioned dimensionedMyTensorN; typedef dimensioned dimensionedMyDiagTensorN; typedef dimensioned dimensionedMySphericalTensorN; // typedef dimensioned dimensionedMySymmTensorN; typedef Foam::Field myVectorNField; typedef GeometricField volMyVectorNField; typedef GeometricField volMyTensorNField; typedef GeometricField volMyDiagTensorNField; typedef GeometricField volMySphericalTensorNField;``` Of course this is not complete at the moment. 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

 July 7, 2011, 04:17 #5 Senior Member   Kathrin Kissling Join Date: Mar 2009 Location: Besigheim, Germany Posts: 134 Rep Power: 10 A follow up... Should we just collect this discussion on the extend-project webpage? We could start a group there with special interest on the block solver! Best Kathrin

 July 7, 2011, 10:32 #6 Senior Member   David Gaden Join Date: Apr 2009 Location: Winnipeg, Canada Posts: 436 Rep Power: 15 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.

 July 7, 2011, 11:28 #7 Senior Member   Kathrin Kissling Join Date: Mar 2009 Location: Besigheim, Germany Posts: 134 Rep Power: 10 Hi David, here are the includes Code: ```#ifndef myNThings_H #define myNThings_H #include "scalar.H" #include "VectorN.H" #include "TensorN.H" #include "contiguous.H"``` The namespace is Foam 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

 July 8, 2011, 09:24 #8 Member   Ivor Clifford Join Date: Mar 2009 Location: Switzerland Posts: 91 Rep Power: 10 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 OF-1.6-ext. Regards Ivor

 July 8, 2011, 12:06 #9 Senior Member   David Gaden Join Date: Apr 2009 Location: Winnipeg, Canada Posts: 436 Rep Power: 15 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 > blockM’ has initializer but incomplete type".` I haven't tried using the neutronDiffusion framework with this... I'll look into this now. 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::eqOp(V1&, const V2&, EqOp) [with V1 = Foam::VectorSpace, double, 1369>, V2 = Foam::VectorSpace, double, 1369>, EqOp = Foam::eqOp, int N = 1369, int I = 193]’: .../src/OpenFOAM/lnInclude/VectorSpaceM.H:27: error: template instantiation depth exceeds maximum of 200 (use -ftemplate-depth-NN to increase the maximum) instantiating ‘VectorSpaceOps<1369, 193>::endLoop’``` I was able to get N=12 to work (the compile took a long time, though)... so I'm wondering if the limit may be my compiler.

 July 8, 2011, 12:22 #10 Member   Ivor Clifford Join Date: Mar 2009 Location: Switzerland Posts: 91 Rep Power: 10 37 variables!!! That's way bigger than I've ever tried. You did the right thing trying to instantiate VectorN. You will probably have to do the same with TensorN 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 -ftemplate-depth-XXX 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.

 July 8, 2011, 12:43 #11 Member   Ivor Clifford Join Date: Mar 2009 Location: Switzerland Posts: 91 Rep Power: 10 Just some more food for thought. You might want to think seriously about whether you need the block-coupled 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 block-coupled 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 Runge-Kutta and predictor-corrector type methods instead. That being said, I'm also interested to see how the block-coupled solver handles this one. Please keep us informed.

 July 8, 2011, 14:38 #12 Senior Member   David Gaden Join Date: Apr 2009 Location: Winnipeg, Canada Posts: 436 Rep Power: 15 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 direct-solving the 37x37 reaction matrix on a cell-by-cell basis across the whole mesh 10 to 1000 times per timestep? My problem requires a point-implicit 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 cell-by-cell 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 ODE-type of solution, but the Courrant number requirement (Co < 1) limits the timestep too much, so it would take years to achieve a solution.

 July 8, 2011, 16:51 #13 Member   Ivor Clifford Join Date: Mar 2009 Location: Switzerland Posts: 91 Rep Power: 10 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 point-implicit 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.

 July 9, 2011, 15:11 #14 Senior Member   Ben K Join Date: Feb 2010 Location: Ottawa, Canada Posts: 140 Rep Power: 12 I don't want to hijack this thread, but are there any examples of the block-coupled 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).

 July 9, 2011, 15:21 #15 Member   Ivor Clifford Join Date: Mar 2009 Location: Switzerland Posts: 91 Rep Power: 10 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.

July 14, 2011, 04:58
non-searchable discussion
#16
Member

Roland
Join Date: Mar 2009
Location: Netherlands
Posts: 65
Rep Power: 10
Hi,

Why move away from an open, Google-searchable discussion to a closed, non-searchable discussion? Don't get me wrong, I like the extend-project, but I don't see the benefit of having an additinal forum/discussion-site which isn't even searchable.

regards,
Sylvester

Quote:
 Originally Posted by kathrin_kissling A follow up... Should we just collect this discussion on the extend-project webpage? We could start a group there with special interest on the block solver! Best Kathrin

 July 14, 2011, 14:08 #17 Senior Member   David Gaden Join Date: Apr 2009 Location: Winnipeg, Canada Posts: 436 Rep Power: 15 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) - Non-searchable (alleged) I'm not sure why the portal isn't searchable. It should be possible to do, given that other database / php-interface websites are (e.g. wikipedia). If the Portal wants to displace the cfd-online 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 Extend-related 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.

July 15, 2011, 07:18
#18
Assistant Moderator

Bernhard Gschaider
Join Date: Mar 2009
Posts: 4,017
Rep Power: 43
Quote:
 Originally Posted by marupio 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) - Non-searchable (alleged) I'm not sure why the portal isn't searchable. It should be possible to do, given that other database / php-interface websites are (e.g. wikipedia). If the Portal wants to displace the cfd-online 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 Extend-related 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.
AFAIK the extend portal has no intention to replace this forum. The forums there are for groups with special interests (I guess only a small fraction of the OF-population is interested in the organisation of the next Austrian User Group Meeting and there is no need to polute the forum with exchanges about an ongoing development that is only of interest for the developers because the code is not yet available/existing - but that kind of stuff would be usually be discussed in a mailing-list anyway). But as soon as the stuff is existing this Forum is the first choice.

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)

 July 15, 2011, 10:44 #19 Senior Member   David Gaden Join Date: Apr 2009 Location: Winnipeg, Canada Posts: 436 Rep Power: 15 Fair enough. The "Community-driven Collaboration" section of the portal looks very functional too.

 July 15, 2011, 11:14 #20 Senior Member   David Gaden Join Date: Apr 2009 Location: Winnipeg, Canada Posts: 436 Rep Power: 15 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 "Community-driven Collaborative Projects" link on the side bar under community. Then choose "Block Solver in OpenFOAM". Files are available from a link on the top.

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post nakul OpenFOAM 7 January 3, 2017 10:25 Peter_600 OpenFOAM 4 August 2, 2014 09:52 nakul OpenFOAM Announcements from Other Sources 18 February 19, 2013 08:48 nishant_hull OpenFOAM Programming & Development 0 August 25, 2009 12:48 bearcat CFX 6 April 28, 2008 14:08

All times are GMT -4. The time now is 19:21.