CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM (https://www.cfd-online.com/Forums/openfoam/)
-   -   blockCoupled solver (https://www.cfd-online.com/Forums/openfoam/90312-blockcoupled-solver.html)

marupio July 6, 2011 18:34

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?

cliffoi July 6, 2011 22:44

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

cliffoi July 6, 2011 22:48

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.

kathrin_kissling July 7, 2011 04:33

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 July 7, 2011 05:17

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

marupio July 7, 2011 11:32

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.

kathrin_kissling July 7, 2011 12:28

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

cliffoi July 8, 2011 10:24

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

marupio July 8, 2011 13:06

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.

cliffoi July 8, 2011 13:22

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 July 8, 2011 13:43

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.

marupio July 8, 2011 15:38

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.

cliffoi July 8, 2011 17:51

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.

benk July 9, 2011 16:11

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).

cliffoi July 9, 2011 16:21

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.

sylvester July 14, 2011 05:58

non-searchable discussion
 
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 (Post 315078)
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


marupio July 14, 2011 15:08

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.

gschaider July 15, 2011 08:18

Quote:

Originally Posted by marupio (Post 316168)
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)

marupio July 15, 2011 11:44

Fair enough. The "Community-driven Collaboration" section of the portal looks very functional too.

marupio July 15, 2011 12:14

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.

holger_marschall July 16, 2011 14:25

Hi Sylvester,
Hi David,

I am the administrator of the Extend Community Portal.

You are indeed making a valid point. However, my view is a bit more differentiated: It is my strong belief that development must be as transparent as possible AND that freedom must be provided to developers to the greatest extend possible. The ultimate goal of providing infrastructure (a content management system in this case: forum, project management, etc.) for the project is to stimulate developments. This **must** be provided both free (of costs) and free (in the sense of freedom).

This said, there are surely many possibilities to approach this goal. The Extend Community Portal provides some; the members are free to choose among different options:
  1. Groups
    There is the option to create Special Interest Groups, Local User Groups and Working Groups. The third 'type', a working group would be about working on a specific issue using OpenFOAM technology. A group forum (that is automatically created) is open and transparent. All forums are searchable.
  2. Project
    This is additional functionality providing the possbility to manage a project (file sharing, task and time management, etc.). A project can be created as closed or open project. Eventually that is up to the founder.

Sylvester, I think, it's really fair to say that with code being in a very early stage (definitely far from a release), people should be allowed to organize their work on their code using free and open accessible infrastructure as provided by Extend. If this is organised around open or closed projects is up to the users and would happen anyway. I would surely like to see that developers **always** create a 'working group' when starting a project on the portal, where they announce that they have started to work on a specific topic, so people can ask questions or ask how to join. Then, the next step would be to do a wall post on the corresponding SIG site being close to the working area. In the past, one would not even been aware of anything going on, because people are sharing e-mails ;). So the portal is putting together many functions to get informations organized **at place** and to provide the possibility to contact the members of a group that are currently active in the respective field.

David, thanks for your input. You are completely right. The CMS this site is build around is a modular one (diffent layout per functionality). One surely must take some time to get into it and I do appreciate any **concrete** hint how to improve structure. For this, please report your suggestions to the bug report system (Mantis) of Extend. It's like the documentation... let's improve it in a collaborative manner! :)

best regards,
Holger


All times are GMT -4. The time now is 16:47.