[deal.II] Axisymmetric problems

[Andreas Kyritsakis]

Dear all,

I have been using Deal.II for quite a time, as the basis of the multi-scale 
multi-physics code FEMOCS <https://github.com/veskem/femocs> we have been 
developing in our groups. Until now our code has been strictly 3D, but now 
we need to solve some problems with rotational symmetry and we want to 
adapt the code for 2D axisymmetric problems in order to save computational 
power.

I have been trying to figure out how to set up such a system, but to be 
honest I am completely lost. I understand that I have to use a different 
mapping than the default one when definig the FEValues object (probably a 
mapping that is associated with a CylindricalManifold, but I cannot figure 
out how to do it, especially how to associate a certain manifold to a 
mapping. 

In general, I think the webpage is really missing a tutorial or an example 
where a simple 2D axisymmetric problem is solved (e.g. just the 
laplace/poisson equation in such a domain). Such problems are very common 
and a lot of people would benefit from such a tutorial. 

Can anyone help in this? 

Thank you very much,
Andreas Kyritsakis

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/c268eea0-19cd-4c3a-93aa-6e866fe2de6d%40googlegroups.com.

Re: [deal.II] Associating a MappingManifold with a Manifold

[Andreas Kyritsakis]

Dear Luca,

Thank you very much for your help. I tried to use your recommendation, but 
it does not seem to work correctly. My code is:
  
CylindricalManifold<2> cylman(0);
triangulation.set_manifold(numbers::flat_manifold_id, cylman);

QGauss<2> quadrature_formula(fe.degree + 1);
MappingManifold<2> mapping;
FEValues<2> fe_values(mapping, fe, quadrature_formula,  update_values | 
update_gradients | update_JxW_values);

The fe_values object seems to have nonsense inside. All the shape_gradients 
are nans and the JxW values are zeros. I also tried to check whether the 
internaldata of the MappingManifold are correct, but they seem to not be 
initialized. I printed

auto  = dynamic_cast< MappingManifold<2>::InternalData 
&>(*mapping.get_data(update_values | update_gradients | update_JxW_values, 
quadrature_formula));
std::cout << "MappingManifold manifold location:" << data.manifold << 
std::endl;
for (auto vol : data.volume_elements)
  std::cout << "volume elements :" << vol << std::endl;

The data.manifold seems to be a NULL while the data.volume_elements is 
empty. I guess I am not initializing something correctly, but I can't find 
what. 

Any help would be very much appreciated.

Best regards,
Andreas



On Tuesday, November 12, 2019 at 8:04:04 PM UTC+2, luca.heltai wrote:
>
> You don’t need to do anything special. Just create a MappingManifold, and 
> pass it where required. The manifold that is used is the one associated to 
> the triangulation. You don’t need to associate it to the Mapping, since the 
> mapping is used on *all* objects of the triangulation, and it will know how 
> to construct the data from the underlying manifold objects. 
>
> Best, 
> Luca. 
>
> > On 12 Nov 2019, at 15:37, Andreas Kyritsakis  > wrote: 
> > 
> > 
> > Hi, 
> > 
> > I want to use a real-to-reference-cell mapping which corresponds to 
> cylindrical symmetry in 2D. In other words, the integration jacobian and 
> the shape function gradients should contain the scale factor of the 
> cylindrical coordinates. I tried to do it with the following code: 
> > 
> >   QGauss<2> quadrature_formula(fe.degree + 1); 
> >   CylindricalManifold<2> cyl_man(0); 
> >   MappingManifold<2> mapping; 
> >   auto  = dynamic_cast< MappingManifold<2>::InternalData 
> &>(*mapping.get_data(update_values | update_gradients | update_JxW_values, 
> quadrature_formula)); 
> >   data.manifold = _man; 
> > 
> > but it does not seem to work. I am afraid that this is not the correct 
> way to associate the Manifold object  with the Mapping object, but I can't 
> figure out how it should be done. 
> > 
> > Can anyone help? 
> > 
> > Best regards, 
> > Andreas 
> > 
> > -- 
> > The deal.II project is located at http://www.dealii.org/ 
> > For mailing list/forum options, see 
> https://groups.google.com/d/forum/dealii?hl=en 
> > --- 
> > You received this message because you are subscribed to the Google 
> Groups "deal.II User Group" group. 
> > To unsubscribe from this group and stop receiving emails from it, send 
> an email to dea...@googlegroups.com . 
> > To view this discussion on the web visit 
> https://groups.google.com/d/msgid/dealii/9d0b3e81-6559-4dc9-bc89-0f4a74ec4921%40googlegroups.com.
>  
>
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/da45ae15-fd2b-4002-aacf-70f3e569bcc2%40googlegroups.com.

Re: [deal.II] Axisymmetric problems

[Andreas Kyritsakis]

Wolfgang,

> How does GetDP do that? Is it a global flag that is set once? 

In GetDP the user defines in the input script an object called "Jacobian", 
which can be different for each domain region. When setting up an integral 
term in the weak formulation then, you have to specify which jacobian 
method will be used to evaluate the integrals. In analogy to Deal.II 
terminology, you define a certain mapping or FEValues object (I am not 
exactly sure which class calculates the jacobian in DealII), which applies 
to different parts of the triangulation (e.g. separated by material_id). 
Then whenever evaluating integrals in a certain part of the triangulation, 
the corresponding mapping/FEValue object would invoke the corresponding 
appropriate JxW that corresponds to that triangulartion object.

> What if one wanted to solve a coupled problem where one PDE is 
axisymmetric but another 
> one is fully 3d? 

In that case the domains/meshes should be different right? One domain is 2D 
the other is 3D. In that case one has just to define the desired jacobian 
to the domain that is axisymmetric, and declare that this jacobian is to be 
used when writing the corresponding weak form terms.



On Thursday, November 14, 2019 at 4:32:37 PM UTC+2, Wolfgang Bangerth wrote:
>
>
> > It is not cumbersome at all to implement it in the assembly function. I 
> > implemented it, and it is just one extra line of code in the assembly 
> > function. The big advantage I find in somehow including that 2πr factor 
> in the 
> > JxW(), applies when having a big code that solves several PDEs on the 
> same 
> > domain, using 3D, 2D or axisymmetric geometries in general. In that 
> case, it 
> > would be advantageous not to need to change every single assembly 
> function 
> > depending on what kind of geometry you have, just use a different 
> FEValues and 
> > implicitly include that factor on the JxW(). This would be in line with 
> the 
> > generic programming approach of Deal.II for dimension-independent 
> > implementations, using the  templates. 
>
> How does GetDP do that? Is it a global flag that is set once? What if one 
> wanted to solve a coupled problem where one PDE is axisymmetric but 
> another 
> one is fully 3d? 
>
>
> > I think it was my big misunderstanding. I misinterpreted some points in 
> this 
> > <
> https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fdealii.org%2Fdeveloper%2Fdoxygen%2Fdeal.II%2FclassMappingManifold_1_1InternalData.html=02%7C01%7CWolfgang.Bangerth%40colostate.edu%7C553194e059364721c87308d768ec7f3e%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637093237276121270=Eke2ZtsNex86tnv%2FVM2gjgYiOYLEFrsklc5GSRnd6Tc%3D=0>
>  
>
> > documentation page, but it was completely my mistake and I think the 
> page is 
> > well-written. I could not understand how to connect the manifold to the 
> > Mapping, but now it is clear to me. Thank you very much for clarifying. 
> > 
> >  > If you wanted to, that would 
> >  > make for a nice pull request to get used to contributing to the 
> library! ;-) ) 
> > 
> > Thank you for the invitation. I will try to see whether and what I can 
> contribute. 
>
> OK, thank you in advance! We're looking forward to it ;-) 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/e67c9539-25b4-4800-b619-a32b9d48d9f8%40googlegroups.com.

[deal.II] Associating a MappingManifold with a Manifold

[Andreas Kyritsakis]

Hi,

I want to use a real-to-reference-cell mapping which corresponds to 
cylindrical symmetry in 2D. In other words, the integration jacobian and 
the shape function gradients should contain the scale factor of the 
cylindrical coordinates. I tried to do it with the following code:

  QGauss<2> quadrature_formula(fe.degree + 1);
  CylindricalManifold<2> cyl_man(0);
  MappingManifold<2> mapping;
  auto  = dynamic_cast< MappingManifold<2>::InternalData 
&>(*mapping.get_data(update_values | update_gradients | update_JxW_values, 
quadrature_formula));
  data.manifold = _man;

but it does not seem to work. I am afraid that this is not the correct way 
to associate the Manifold object  with the Mapping object, but I can't 
figure out how it should be done.

Can anyone help?

Best regards,
Andreas

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/9d0b3e81-6559-4dc9-bc89-0f4a74ec4921%40googlegroups.com.

Re: [deal.II] Axisymmetric problems

[Andreas Kyritsakis]

Wolfgang, 

Thank you very much for your prompt response.

You are right that one could just leave the mapping as a simple 2D and 
change the weak formulation. However I think that it is possible to solve 
an axisymmetric problem without touching the weak formulation. I have been 
working for quite some time with GetDP, and there that is the standard way. 
For axisymmetric problems one just needs to change the "Jacobian" which in 
GetDP stands for the mapping between real cells and the reference/unit 
cell. I think the two approaches in the end boil down in the same equations 
and system assembly. In one case the coordinate transform scaling factors 
are considered as a part of the PDE and need to be explicitly inserted in 
the weak form, while in the other they are considered as part of the cell 
mapping and would appear implicitly in the jacobian of the transfrom (the 
JxW factors in the Deal.II framework). 

>From a programming point of view, I think that the latter approach would 
produce more generic and reusable code, that is why I first tried to pursue 
it. I thought that this could be done in Deal.II, similarly to GetDP, by 
using an appropriate Mapping object. I guess such a Mapping class is not 
implemented though. In case I am wrong and it is (maybe a MappingManifold 
associated with a CylindricalManifold), would it be possible to help on how 
to implement it? 

In general, how can I associate a MappingManifold object with a certain 
Manifold object? In the documentation, such an association is mentioned 
several times, but it is not specified how it can be implemented.

Cheers,
Andreas


On Thursday, November 7, 2019 at 7:40:37 PM UTC+2, Wolfgang Bangerth wrote:
>
>
> Andreas, 
>
> > I have been using Deal.II for quite a time, as the basis of the 
> multi-scale 
> > multi-physics code FEMOCS 
> > <
> https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fveskem%2Ffemocs=02%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cd7d8d7ed6a644b803f8908d7639a4fb1%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C1%7C637087386726073424=oXStQieIrHcNMPw4sO%2FfnCF%2B8mxdRgGXSrprxUrMo0w%3D=0>
>  
>
> > we have been developing in our groups. Until now our code has been 
> strictly 
> > 3D, but now we need to solve some problems with rotational symmetry and 
> we 
> > want to adapt the code for 2D axisymmetric problems in order to save 
> > computational power. 
> > 
> > I have been trying to figure out how to set up such a system, but to be 
> honest 
> > I am completely lost. I understand that I have to use a different 
> mapping than 
> > the default one when definig the FEValues object (probably a mapping 
> that is 
> > associated with a CylindricalManifold, but I cannot figure out how to do 
> it, 
> > especially how to associate a certain manifold to a mapping. 
>
> No, if you use cylindrical coordinates, then you work in r-z space and 
> your 
> domain is probably just a rectangle. But your PDE changes. 
>
> As for resources: Yours is a rather frequently asked question. Rather than 
> trying to repeat what I've been writing on many other occasions, have you 
> tried to search through the mailing list archives to find some starting 
> point? 
>
>
> > In general, I think the webpage is really missing a tutorial or an 
> example 
> > where a simple 2D axisymmetric problem is solved (e.g. just the 
> > laplace/poisson equation in such a domain). Such problems are very 
> common and 
> > a lot of people would benefit from such a tutorial. 
>
> Yes, we should eventually write such a tutorial... 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/b3127fda-f4b5-4d01-b206-c2077b50d431%40googlegroups.com.

[deal.II] Bug Report in GridTools::get_active_neighbors()

[Andreas Kyritsakis]

Hello,

This is to report what I think is a bug in 
GridTools::get_active_neighbors().

The following code should compile but it does not.

#include 
#include 
#include 
#include 
using namespace dealii;

template
class Test{
public: 
Test(Triangulation *tria) : triangulation(triangulation), 
dof_handler(*tria){}

void run(){
typename DoFHandler::active_cell_iterator cell;
  
for (cell = dof_handler.begin_active(); cell != dof_handler.end(); 
++cell){
std::vector::active_cell_iterator> 
neighbors;
GridTools::get_active_neighbors(cell, neighbors); 
}
}
private:
Triangulation* triangulation;
DoFHandler dof_handler;
};

int main(){
Triangulation<2> triangulation;
GridGenerator::hyper_cube(triangulation, -1, 1);
Test<2> test();
test.run();
}

I think somehow the compiler cannot resolve the generic type 
MeshType::active_cell_iterator into the specific one 
DoFHandler::active_cell_iterator, although from what I understand in 
the documentation it should. Probably some template specification missing 
somewhere, but I can't find it.

In my code I fixed it easily by copying the function directly to my own 
files and specifying appropriately, but I think it would be nice to report 
the bug.

Cheers,
Andreas 

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/27895150-8268-4f88-96e9-634172564a44%40googlegroups.com.

Re: [deal.II] Axisymmetric problems

[Andreas Kyritsakis]

Hi Wolfgang,

Thank you very much. Your answer was very helpful. 

> You're correct -- the mapping classes don't do that. I also think that 
they
> shouldn't: They have a different purpose and at least at the moment I 
don't
> quite see how that could be put into what the mappings are supposed to 
do. The
> place to do it might be where the FEValues class combines quadrature 
weights
> with the determinant of the Jacobian -- though I might still think that 
it's
> not so cumbersome to just keep everything in user code -- it'll be at 
most a
> dozen lines of code for most applications.

It is not cumbersome at all to implement it in the assembly function. I 
implemented it, and it is just one extra line of code in the assembly 
function. The big advantage I find in somehow including that 2πr factor in 
the JxW(), applies when having a big code that solves several PDEs on the 
same domain, using 3D, 2D or axisymmetric geometries in general. In that 
case, it would be advantageous not to need to change every single assembly 
function depending on what kind of geometry you have, just use a different 
FEValues and implicitly include that factor on the JxW(). This would be in 
line with the generic programming approach of Deal.II for 
dimension-independent implementations, using the  templates. 

> Can you point out where the places you mention about the documentation 
are?
> Maybe we can clarify the text there.

I think it was my big misunderstanding. I misinterpreted some points in this 

 
documentation page, but it was completely my mistake and I think the page 
is well-written. I could not understand how to connect the manifold to the 
Mapping, but now it is clear to me. Thank you very much for clarifying.

> If you wanted to, that would
> make for a nice pull request to get used to contributing to the library! 
;-) )

Thank you for the invitation. I will try to see whether and what I can 
contribute.

On Thursday, November 14, 2019 at 6:02:33 AM UTC+2, Wolfgang Bangerth wrote:
>
>
> Andreas, 
>
> > You are right that one could just leave the mapping as a simple 2D and 
> change 
> > the weak formulation. However I think that it is possible to solve an 
> > axisymmetric problem without touching the weak formulation. 
>
> You can, in principle, but it requires some tricks. If you want to use an 
> axisymmetric formulation, then in essence you formulate your shape 
> functions as 
>phi_i(r,phi,z) 
> so that they are only defined on a r/z mesh, with no dependence on phi -- 
> but 
> nevertheless defined in 3d. The point is that you have to have a 2d mesh 
> for 
> the cross section in 3d. It's not a convenient set up. It's easier to 
> recall 
> that if you have 
>
>\int_\Omega   F(r(x,y,z),z)  dx dy dz 
>
> (i.e., with an integrand that formally depends on x,y,z, but in reality 
> only 
> on r,z) can be expressed as 
>
>\int_a^b \int_0^{2pi} \int_0^R   F(r,z)  r dr dphi dz 
>=   2pi \int_a^b \int_0^R   F(r,z)  r dr dz 
>
> In other words, the only change to the bilinear form is the addition of a 
> factor of 2*pi*r to all integrals that might appear in the formulation -- 
> not 
> a bothersome task. 
>
>
> > I have been 
> > working for quite some time with GetDP, and there that is the standard 
> way. 
>
> I'm curious how that works there. 
>
>
> > For axisymmetric problems one just needs to change the "Jacobian" which 
> in 
> > GetDP stands for the mapping between real cells and the reference/unit 
> cell. 
>
> So maybe the only difference there is the factor of 2*pi*r to |det J| 
> (which 
> is part of the JxW factor that appears in every deal.II integral). 
>
>
> > I 
> > think the two approaches in the end boil down in the same equations and 
> system 
> > assembly. In one case the coordinate transform scaling factors are 
> considered 
> > as a part of the PDE and need to be explicitly inserted in the weak 
> form, 
> > while in the other they are considered as part of the cell mapping and 
> would 
> > appear implicitly in the jacobian of the transfrom (the JxW factors in 
> the 
> > Deal.II framework). 
>
> Yes, that seems reasonable. 
>
>
> >  From a programming point of view, I think that the latter approach 
> would 
> > produce more generic and reusable code, that is why I first tried to 
> pursue 
> > it. I thought that this could be done in Deal.II, similarly to GetDP, by 
> using 
> > an appropriate Mapping object. I guess such a Mapping class is not 
> implemented 
> > though. In case I am wrong and it is (maybe a MappingManifold associated 
> with 
> > a CylindricalManifold), would it be possible to help on how to implement 
> it? 
>
> You're correct -- the mapping classes don't do that. I also think that 
> they 
> shouldn't: They have a different purpose and at least at the moment I 
> don't 
> quite see how that could be put into what the mappings are supposed to do. 
> The 
> place to 

[deal.II] problem measure() for quads

[Andreas Kyritsakis]

Hi,

I have a problem with using the measure() function for quads in 3D. I am 
currently using version 9.1.1 and the function seems to check whether the 
face is planar and if not, it returns nan. The problem is that the 
tolerance for considering it planar is extremely low (1.e-24), which causes 
problems. If for example the mesh is read in from a standard .msh file, 
written by the GridOut class, we cannot expect to get such a precision in 
the vertices. I saw that in the 9.2 version this is partially fixed by 
handling the non-planar case by numerical intergration, but the tolerance 
has remained the same extremely low, which might result in unnecessary 
computational effort. 

I wonder whether the tolerance can be increased to something more 
reasonable (e.g. 1.e-10) or there is something I am missing that suggests 
such a low tolerance.

Best regards,
Andreas Kyritsakis

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/b52ebd60-f264-4512-92f5-144f1e3d5f5a%40googlegroups.com.

Re: [deal.II] problem measure() for quads

[Andreas Kyritsakis]

Thank you very much Wolfgang. Yes I would like to propose a patch. How is 
it done? I am not familiar with the process.

On Thursday, May 28, 2020 at 7:31:04 PM UTC+3, Wolfgang Bangerth wrote:
>
> On 5/28/20 3:08 AM, Andreas Kyritsakis wrote: 
> > 
> > I have a problem with using the measure() function for quads in 3D. I am 
> > currently using version 9.1.1 and the function seems to check whether 
> the face 
> > is planar and if not, it returns nan. The problem is that the tolerance 
> for 
> > considering it planar is extremely low (1.e-24), which causes problems. 
> If for 
> > example the mesh is read in from a standard .msh file, written by the 
> GridOut 
> > class, we cannot expect to get such a precision in the vertices. I saw 
> that in 
> > the 9.2 version this is partially fixed by handling the non-planar case 
> by 
> > numerical intergration, but the tolerance has remained the same 
> extremely low, 
> > which might result in unnecessary computational effort. 
> > 
> > I wonder whether the tolerance can be increased to something more 
> reasonable 
> > (e.g. 1.e-10) or there is something I am missing that suggests such a 
> low 
> > tolerance. 
>
> It's a good question. I checked that the expression in question is 
> dimensionless and scaled correctly so that the mesh size h does not appear 
> in 
> the size of the expression. I think it would be quite ok if you changed 
> the 
> tolerance to 1e-10. 
>
> Want to propose a patch for this? 
>
> Best 
>   W. 
>
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/231d1ecd-03a2-4cce-9590-ed1470fd5278%40googlegroups.com.

[deal.II] Re: Resolving with different RHS

[Andreas Kyritsakis]

Thanks a lot!!

On Tuesday, June 2, 2020 at 1:21:37 PM UTC+3, Simon Sticko wrote:
>
> Hi,
> a better option than computing the inverse is to factorize the matrix. 
> This can be done using the SparseDirectUMFPACK solver:
>
>
> https://www.dealii.org/current/doxygen/deal.II/classSparseDirectUMFPACK.html
>
> You might want to take a look at step 22, which uses this solver:
>
> https://www.dealii.org/current/doxygen/deal.II/step_22.html
>
> Best,
> Simon
>
> On Tuesday, June 2, 2020 at 12:04:38 PM UTC+2, Andreas Kyritsakis wrote:
>>
>> Dear all,
>>
>> I have a problem where a Poisson-like equation has to be solved again and 
>> again in many time steps. At each timestep the LHS remains the same while 
>> the RHS changes (slightly). My current implementation is to use the 
>> standard SolverCG, passing the old solution as initial, which already 
>> reduces the number of CG steps required. Yet, I wonder whether my approach 
>> is naive and there is a faster way to implement this, taking better 
>> advantage of the fact that the LHS stays always the same. I initially 
>> thought about calculating the inverse matrix once and just doing matrix 
>> multiplication, but since the mass matrix is a huge sparse matrix, only 
>> storing the inverse (which in general is not sparse) would require huge 
>> memory. I also thought about the LinearOperator concepts, but if I 
>> understood correctly, they just implement a nice wrapper to call a solver 
>> each time. Am I missing something?
>>
>> Cheers,
>> Andreas
>>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/072b1775-eaae-486e-a874-2fba0ca4e36d%40googlegroups.com.

[deal.II] Resolving with different RHS

[Andreas Kyritsakis]

Dear all,

I have a problem where a Poisson-like equation has to be solved again and 
again in many time steps. At each timestep the LHS remains the same while 
the RHS changes (slightly). My current implementation is to use the 
standard SolverCG, passing the old solution as initial, which already 
reduces the number of CG steps required. Yet, I wonder whether my approach 
is naive and there is a faster way to implement this, taking better 
advantage of the fact that the LHS stays always the same. I initially 
thought about calculating the inverse matrix once and just doing matrix 
multiplication, but since the mass matrix is a huge sparse matrix, only 
storing the inverse (which in general is not sparse) would require huge 
memory. I also thought about the LinearOperator concepts, but if I 
understood correctly, they just implement a nice wrapper to call a solver 
each time. Am I missing something?

Cheers,
Andreas

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/acd95fd2-612b-4b06-97a2-a72e3584a508%40googlegroups.com.

Re: [deal.II] Resolving with different RHS

[Andreas Kyritsakis]

Thank you both!

On Tuesday, June 2, 2020 at 10:12:43 PM UTC+3, Jean-Paul Pelteret wrote:
>
> Hi Andreas,
>
> I concur with Simon. Additionally, if you’re using the Trilinos linear 
> algebra classes, then the TrilinosWrappers::SolverDirect solver will do the 
> same
>
> https://dealii.org/developer/doxygen/deal.II/classTrilinosWrappers_1_1SolverDirect.html#a6c0fecc921bc93f2c005e265cfbdd882
> Its hard to tell (from the documentation) what 
> PETScWrappers::SparseDirectMUMPS does, but if it doesn’t retain the 
> factorisation then I’m sure some light refactoring of the class could 
> address that issue.
>
> Best,
> Jean-Paul
>
> On 02 Jun 2020, at 12:53, Andreas Kyritsakis  > wrote:
>
> Thanks a lot!!
>
> On Tuesday, June 2, 2020 at 1:21:37 PM UTC+3, Simon Sticko wrote:
>>
>> Hi,
>> a better option than computing the inverse is to factorize the matrix. 
>> This can be done using the SparseDirectUMFPACK solver:
>>
>>
>> https://www.dealii.org/current/doxygen/deal.II/classSparseDirectUMFPACK.html
>>
>> You might want to take a look at step 22, which uses this solver:
>>
>> https://www.dealii.org/current/doxygen/deal.II/step_22.html
>>
>> Best,
>> Simon
>>
>> On Tuesday, June 2, 2020 at 12:04:38 PM UTC+2, Andreas Kyritsakis wrote:
>>>
>>> Dear all,
>>>
>>> I have a problem where a Poisson-like equation has to be solved again 
>>> and again in many time steps. At each timestep the LHS remains the same 
>>> while the RHS changes (slightly). My current implementation is to use the 
>>> standard SolverCG, passing the old solution as initial, which already 
>>> reduces the number of CG steps required. Yet, I wonder whether my approach 
>>> is naive and there is a faster way to implement this, taking better 
>>> advantage of the fact that the LHS stays always the same. I initially 
>>> thought about calculating the inverse matrix once and just doing matrix 
>>> multiplication, but since the mass matrix is a huge sparse matrix, only 
>>> storing the inverse (which in general is not sparse) would require huge 
>>> memory. I also thought about the LinearOperator concepts, but if I 
>>> understood correctly, they just implement a nice wrapper to call a solver 
>>> each time. Am I missing something?
>>>
>>> Cheers,
>>> Andreas
>>>
>>
> -- 
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see 
> https://groups.google.com/d/forum/dealii?hl=en
> --- 
> You received this message because you are subscribed to the Google Groups 
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an 
> email to dea...@googlegroups.com .
> To view this discussion on the web visit 
> https://groups.google.com/d/msgid/dealii/072b1775-eaae-486e-a874-2fba0ca4e36d%40googlegroups.com
>  
> <https://groups.google.com/d/msgid/dealii/072b1775-eaae-486e-a874-2fba0ca4e36d%40googlegroups.com?utm_medium=email_source=footer>
> .
>
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/c1d366a1-fdc2-4798-9a10-dd562fd342f0%40googlegroups.com.