CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM

blockCoupled solver

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

Like Tree1Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   July 6, 2011, 17:34
Default blockCoupled solver
  #1
Senior Member
 
David Gaden
Join Date: Apr 2009
Location: Winnipeg, Canada
Posts: 397
Rep Power: 12
marupio is on a distinguished road
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?
marupio is offline   Reply With Quote

Old   July 6, 2011, 21:44
Default
  #2
Member
 
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 8
cliffoi is on a distinguished road
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>
class mySystemOfEquations
{
    typedef VectorN<cmpt, length> vectorType;
...
    blockSolverPerformance<vectorType> 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<scalar, 2> eqns(...);
        eqns.solve();
    case(4):
        mySystemOfEquations<scalar, 4> 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.
cliffoi is offline   Reply With Quote

Old   July 6, 2011, 21:48
Default
  #3
Member
 
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 8
cliffoi is on a distinguished road
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.
cliffoi is offline   Reply With Quote

Old   July 7, 2011, 03:33
Default
  #4
Senior Member
 
Kathrin Kissling
Join Date: Mar 2009
Location: Besigheim, Germany
Posts: 134
Rep Power: 8
kathrin_kissling is on a distinguished road
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<scalar, N> myVectorN;
    typedef TensorN<scalar, N> myTensorN;
    typedef DiagTensorN<scalar, N> myDiagTensorN;
    typedef SphericalTensorN<scalar, N> mySphericalTensorN;
//     typedef symmTensorN<scalar, N> mySymmTensorN;

    typedef dimensioned<myTensorN> dimensionedMyTensorN;
    typedef dimensioned<myDiagTensorN> dimensionedMyDiagTensorN;
    typedef dimensioned<mySphericalTensorN> dimensionedMySphericalTensorN;
//     typedef dimensioned<mySymmTensorN> dimensionedMySymmTensorN;
    
    typedef Foam::Field<myVectorN> myVectorNField;
    
    typedef GeometricField<myVectorN, fvPatchField, volMesh> volMyVectorNField;
    typedef GeometricField<myTensorN, fvPatchField, volMesh> volMyTensorNField;
    typedef GeometricField<myDiagTensorN, fvPatchField, volMesh> volMyDiagTensorNField;
    typedef GeometricField<mySphericalTensorN, fvPatchField, volMesh> 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
kathrin_kissling is offline   Reply With Quote

Old   July 7, 2011, 04:17
Default
  #5
Senior Member
 
Kathrin Kissling
Join Date: Mar 2009
Location: Besigheim, Germany
Posts: 134
Rep Power: 8
kathrin_kissling is on a distinguished road
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
kathrin_kissling is offline   Reply With Quote

Old   July 7, 2011, 10:32
Default
  #6
Senior Member
 
David Gaden
Join Date: Apr 2009
Location: Winnipeg, Canada
Posts: 397
Rep Power: 12
marupio is on a distinguished road
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.
marupio is offline   Reply With Quote

Old   July 7, 2011, 11:28
Default
  #7
Senior Member
 
Kathrin Kissling
Join Date: Mar 2009
Location: Besigheim, Germany
Posts: 134
Rep Power: 8
kathrin_kissling is on a distinguished road
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
kathrin_kissling is offline   Reply With Quote

Old   July 8, 2011, 09:24
Default
  #8
Member
 
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 8
cliffoi is on a distinguished road
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
cliffoi is offline   Reply With Quote

Old   July 8, 2011, 12:06
Default
  #9
Senior Member
 
David Gaden
Join Date: Apr 2009
Location: Winnipeg, Canada
Posts: 397
Rep Power: 12
marupio is on a distinguished road
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 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<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]’:
.../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.
marupio is offline   Reply With Quote

Old   July 8, 2011, 12:22
Default
  #10
Member
 
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 8
cliffoi is on a distinguished road
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 -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.
cliffoi is offline   Reply With Quote

Old   July 8, 2011, 12:43
Default
  #11
Member
 
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 8
cliffoi is on a distinguished road
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.
cliffoi is offline   Reply With Quote

Old   July 8, 2011, 14:38
Default
  #12
Senior Member
 
David Gaden
Join Date: Apr 2009
Location: Winnipeg, Canada
Posts: 397
Rep Power: 12
marupio is on a distinguished road
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.
marupio is offline   Reply With Quote

Old   July 8, 2011, 16:51
Default
  #13
Member
 
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 8
cliffoi is on a distinguished road
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.
cliffoi is offline   Reply With Quote

Old   July 9, 2011, 15:11
Default
  #14
Senior Member
 
Ben K
Join Date: Feb 2010
Location: Ottawa, Canada
Posts: 140
Rep Power: 10
benk is on a distinguished road
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).
benk is offline   Reply With Quote

Old   July 9, 2011, 15:21
Default
  #15
Member
 
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 8
cliffoi is on a distinguished road
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.
cliffoi is offline   Reply With Quote

Old   July 14, 2011, 04:58
Default non-searchable discussion
  #16
Member
 
Roland
Join Date: Mar 2009
Location: Netherlands
Posts: 62
Rep Power: 8
sylvester is on a distinguished road
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 View Post
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
sylvester is offline   Reply With Quote

Old   July 14, 2011, 14:08
Default
  #17
Senior Member
 
David Gaden
Join Date: Apr 2009
Location: Winnipeg, Canada
Posts: 397
Rep Power: 12
marupio is on a distinguished road
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.
marupio is offline   Reply With Quote

Old   July 15, 2011, 07:18
Default
  #18
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,912
Rep Power: 40
gschaider will become famous soon enoughgschaider will become famous soon enough
Quote:
Originally Posted by marupio View Post
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)
gschaider is offline   Reply With Quote

Old   July 15, 2011, 10:44
Default
  #19
Senior Member
 
David Gaden
Join Date: Apr 2009
Location: Winnipeg, Canada
Posts: 397
Rep Power: 12
marupio is on a distinguished road
Fair enough. The "Community-driven Collaboration" section of the portal looks very functional too.
marupio is offline   Reply With Quote

Old   July 15, 2011, 11:14
Default
  #20
Senior Member
 
David Gaden
Join Date: Apr 2009
Location: Winnipeg, Canada
Posts: 397
Rep Power: 12
marupio is on a distinguished road
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.
marupio is offline   Reply With Quote

Reply

Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
thobois class engineTopoChangerMesh error Peter_600 OpenFOAM 4 August 2, 2014 09:52
A New Solver for Supersonic Combustion nakul OpenFOAM 6 March 30, 2014 03:00
A New Solver for Supersonic Combustion nakul OpenFOAM Announcements from Other Sources 18 February 19, 2013 08:48
Development of a Low mach PISO solver nishant_hull OpenFOAM Programming & Development 0 August 25, 2009 12:48
why the solver reject it? Anyone with experience? bearcat CFX 6 April 28, 2008 14:08


All times are GMT -4. The time now is 23:42.