[petsc-users] multigrid preconditioning and adaptivity

Tue Mar 15 04:43:16 CDT 2016

Hello Matt,

On 7 Mar 2016, at 16:42, Matthew Knepley <knepley at gmail.com<mailto:knepley at gmail.com>> wrote:

Yes, that should be it. It would be nice to have some example that does this if you would
be willing to contribute some version of your code.

Following discussion from last week, I had time to implement MG via DM, see following code,
http://cdash.eng.gla.ac.uk/mofem/_p_c_m_g_set_up_via_approx_orders_8hpp_source.html

Minimal setup for DM and MG looks like that
132<http://cdash.eng.gla.ac.uk/mofem/_p_c_m_g_set_up_via_approx_orders_8hpp.html#a1c1639cf5baf9432f4f78094171a4358> PetscErrorCode DMCreate_MGViaApproxOrders<http://cdash.eng.gla.ac.uk/mofem/_p_c_m_g_set_up_via_approx_orders_8cpp.html#a1c1639cf5baf9432f4f78094171a4358>(DM dm) {
133  PetscErrorCode ierr;
134  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
135  PetscFunctionBegin;
136  if(!dm->data) {
137  dm->data = new DMMGViaApproxOrdersCtx<http://cdash.eng.gla.ac.uk/mofem/struct_d_m_m_g_via_approx_orders_ctx.html>();
138  } else {
139  ((DMCtx<http://cdash.eng.gla.ac.uk/mofem/struct_mo_f_e_m_1_1_d_m_ctx.html>*)(dm->data))->referenceNumber++;
140  }
141  ierr = DMSetOperators_MoFEM<http://cdash.eng.gla.ac.uk/mofem/group__dm.html#ga022055f7270e8fb02e780268a967c24d>(dm); CHKERRQ(ierr);
142  dm->ops->creatematrix = DMCreateMatrix_MGViaApproxOrders<http://cdash.eng.gla.ac.uk/mofem/group__dm.html#gaac13a65e2603262490a3c07cf846bafc>;
143  dm->ops->createglobalvector = DMCreateGlobalVector_MGViaApproxOrders<http://cdash.eng.gla.ac.uk/mofem/_p_c_m_g_set_up_via_approx_orders_8cpp.html#aa2a7f489f973858a921260d1af52d906>;
144  dm->ops->coarsen = DMCoarsen_MGViaApproxOrders<http://cdash.eng.gla.ac.uk/mofem/group__dm.html#ga04beeab0e92e7f209334adb347353090>;
145  dm->ops->createinterpolation = DMCreateInterpolation_MGViaApproxOrders<http://cdash.eng.gla.ac.uk/mofem/_p_c_m_g_set_up_via_approx_orders_8cpp.html#ad3bc29f7e017e757597de9af3f82202b>;
146  ierr = DMKSPSetComputeOperators(dm,ksp_set_operators,NULL); CHKERRQ(ierr);
147  PetscInfo1(dm,"Create DMMGViaApproxOrders reference = %d\n",((DMCtx<http://cdash.eng.gla.ac.uk/mofem/struct_mo_f_e_m_1_1_d_m_ctx.html>*)(dm->data))->referenceNumber);
148  PetscFunctionReturn(0);
149 }

It is not yet perfect, I need to add additional interface functions to improve functionality. If this in any way would be useful for you, pleas feel free to use it. I building now two examples, f.e. http://cdash.eng.gla.ac.uk/mofem/elasticity_8cpp_source.html, see

Register/Create and SetUp DM;
303  DMType dm_name = "ELASTIC_PROB";
304  ierr = DMRegister_MGViaApproxOrders<http://cdash.eng.gla.ac.uk/mofem/group__dm.html#ga1f0d8d905f5e398772b076d9582f2bd8>(dm_name); CHKERRQ(ierr);
305  // ierr = DMRegister_MoFEM(dm_name); CHKERRQ(ierr);
306
307  //craete dm instance
308  DM dm;
309  ierr = DMCreate(PETSC_COMM_WORLD,&dm);CHKERRQ(ierr);
310  ierr = DMSetType(dm,dm_name);CHKERRQ(ierr);
311  //set dm datastruture whict created mofem datastructures
312  ierr = DMMoFEMCreateMoFEM<http://cdash.eng.gla.ac.uk/mofem/group__dm.html#ga95264c2ed2c7295fa52072508b60c5fb>(dm,&m_field,dm_name,bit_level0); CHKERRQ(ierr);
313  ierr = DMSetFromOptions(dm); CHKERRQ(ierr);
314  ierr = DMMoFEMSetIsPartitioned<http://cdash.eng.gla.ac.uk/mofem/group__dm.html#gafaa90d67e85cf3c75096fa17a4575901>(dm,is_partitioned); CHKERRQ(ierr);
315  //add elements to dm
316  ierr = DMMoFEMAddElement<http://cdash.eng.gla.ac.uk/mofem/group__dm.html#gaaec6e1adbf45c2c7fd67548674c1bce5>(dm,"ELASTIC"); CHKERRQ(ierr);
317  ierr = DMMoFEMAddElement<http://cdash.eng.gla.ac.uk/mofem/group__dm.html#gaaec6e1adbf45c2c7fd67548674c1bce5>(dm,"BODY_FORCE"); CHKERRQ(ierr);
318  ierr = DMMoFEMAddElement<http://cdash.eng.gla.ac.uk/mofem/group__dm.html#gaaec6e1adbf45c2c7fd67548674c1bce5>(dm,"FLUID_PRESSURE_FE"); CHKERRQ(ierr);
319  ierr = DMMoFEMAddElement<http://cdash.eng.gla.ac.uk/mofem/group__dm.html#gaaec6e1adbf45c2c7fd67548674c1bce5>(dm,"FORCE_FE"); CHKERRQ(ierr);
320  ierr = DMMoFEMAddElement<http://cdash.eng.gla.ac.uk/mofem/group__dm.html#gaaec6e1adbf45c2c7fd67548674c1bce5>(dm,"PRESSURE_FE"); CHKERRQ(ierr);
321  ierr = DMSetUp(dm); CHKERRQ(ierr);

Matrices and vector via DM,
326  Vec F,D,D0;
327  ierr = DMCreateGlobalVector(dm,&F); CHKERRQ(ierr);
328  ierr = VecDuplicate(F,&D); CHKERRQ(ierr);
329  ierr = VecDuplicate(F,&D0); CHKERRQ(ierr);
330  Mat Aij;
331  ierr = DMCreateMatrix(dm,&Aij); CHKERRQ(ierr);

SetUP solver,
421  KSP solver;
422  ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(ierr);
423  ierr = KSPSetDM(solver,dm); CHKERRQ(ierr);
424  ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
425  ierr = KSPSetOperators(solver,Aij,Aij); CHKERRQ(ierr);
426  {
427  //from PETSc example ex42.c
428  PetscBool same = PETSC_FALSE;
429  PC pc;
430  ierr = KSPGetPC(solver,&pc); CHKERRQ(ierr);
431  PetscObjectTypeCompare((PetscObject)pc,PCMG,&same);
432  if (same) {
433  PCMGSetUpViaApproxOrdersCtx<http://cdash.eng.gla.ac.uk/mofem/struct_p_c_m_g_set_up_via_approx_orders_ctx.html> pc_ctx(dm,Aij);
434  ierr = PCMGSetUpViaApproxOrders<http://cdash.eng.gla.ac.uk/mofem/_p_c_m_g_set_up_via_approx_orders_8cpp.html#a86375678c4c815cecd1931d7f51e40c3>(pc,&pc_ctx); CHKERRQ(ierr);
435  } else {
436  // Operators are already set, do not use DM for doing that
437  ierr = KSPSetDMActive(solver,PETSC_FALSE); CHKERRQ(ierr);
438  }
439  }
440  ierr = KSPSetInitialGuessKnoll(solver,PETSC_FALSE); CHKERRQ(ierr);
441  ierr = KSPSetInitialGuessNonzero(solver,PETSC_TRUE); CHKERRQ(ierr);
442  ierr = KSPSetUp(solver); CHKERRQ(ierr);

Kind regards,
Lukasz

On 7 Mar 2016, at 17:41, Lukasz Kaczmarczyk <Lukasz.Kaczmarczyk at glasgow.ac.uk<mailto:Lukasz.Kaczmarczyk at glasgow.ac.uk>> wrote:

On 7 Mar 2016, at 16:42, Matthew Knepley <knepley at gmail.com<mailto:knepley at gmail.com>> wrote:

On Mon, Mar 7, 2016 at 10:28 AM, Lukasz Kaczmarczyk <Lukasz.Kaczmarczyk at glasgow.ac.uk<mailto:Lukasz.Kaczmarczyk at glasgow.ac.uk>> wrote:
Many thanks all for help. I started to implement function for DM.

I understand that minimal implementation is that for the DM i need to have, is to have  DMCoarsen and in each level for all DMs, set operators
DMKSPSetComputeOperators and DMCreateInterpolation. Matrix matrix free P from DMCreateInterpolation have to have operators for mult and mult_traspose. Is that is all?

Yes, that should be it. It would be nice to have some example that does this if you would
be willing to contribute some version of your code.

No problem, I will do that will pleasure.

Kind regards,
Lukasz

On 7 Mar 2016, at 15:55, Mark Adams <mfadams at lbl.gov<mailto:mfadams at lbl.gov>> wrote:

You can just set the coarse grid matrix/operator instead of using Galerkin.  If you have a shell (matrix free) P then you will need to create and set this yourself.  Our Galerkin requires a matrix P.

On Mon, Mar 7, 2016 at 9:32 AM, Lukasz Kaczmarczyk <Lukasz.Kaczmarczyk at glasgow.ac.uk<mailto:Lukasz.Kaczmarczyk at glasgow.ac.uk>> wrote:

> On 7 Mar 2016, at 14:21, Lawrence Mitchell <lawrence.mitchell at imperial.ac.uk<mailto:lawrence.mitchell at imperial.ac.uk>> wrote:
>
> On 07/03/16 14:16, Lukasz Kaczmarczyk wrote:
>>
>>> On 7 Mar 2016, at 13:50, Matthew Knepley <knepley at gmail.com<mailto:knepley at gmail.com>
>>> <mailto:knepley at gmail.com<mailto:knepley at gmail.com>>> wrote:
>>>
>>> On Mon, Mar 7, 2016 at 6:58 AM, Lukasz Kaczmarczyk
>>> <Lukasz.Kaczmarczyk at glasgow.ac.uk<mailto:Lukasz.Kaczmarczyk at glasgow.ac.uk>
>>> <mailto:Lukasz.Kaczmarczyk at glasgow.ac.uk<mailto:Lukasz.Kaczmarczyk at glasgow.ac.uk>>> wrote:
>>>
>>>    Hello,
>>>
>>>    I run multi-grid solver, with adaptivity, works well, however It
>>>    is some space for improving efficiency. I using hierarchical
>>>    approximation basis, for which
>>>    construction of interpolation operators is simple, it is simple
>>>    injection.
>>>
>>>    After each refinement level (increase of order of approximation
>>>    on some element) I rebuild multigrid pre-conditioner with
>>>    additional level. It is a way to add dynamically new levels
>>>    without need of rebuilding whole MG pre-conditioner.
>>>
>>>
>>> That does not currently exist, however it would not be hard to add,
>>> since the MG structure jsut consists of
>>> arrays of pointers.
>>>
>>>
>>>    Looking at execution profile I noticed that 50%-60% of time is
>>>    spent on MatPtAP function during PCSetUP stage.
>>>
>>>
>>> Which means you are using a Galerkin projection to define the coarse
>>> operator. Do you have a direct way of defining
>>> this operator (rediscretization)?
>>
>> Matt,
>>
>>
>> Thanks for swift response. You are right, I using Galerkin projection.
>>
>> Yes, I have a way to get directly coarse operator, it is some sub
>> matrix of whole matrix. I taking advantage here form hierarchical
>> approximation.
>>
>> I could reimplement PCSetUp_MG to set the MG structure directly, but
>> this probably not good approach, since my implementation which will
>> work with current petsc version could be incompatible which  future
>> changes in native MG data structures. The alternative option is to
>> hack MatPtAP itself, and until petsc MG will use this, whatever
>> changes you will make in MG in the future my code will work.
>
> Why not provide a shell DM to the KSP that knows how to compute the
> operators (and how to refine/coarsen and therefore
> restrict/interpolate).  Then there's no need to use Galerkin coarse
> grid operators, and the KSP will just call back to your code to create
> the appropriate matrices.

Hello Lawrence,

Thanks, it is good advice.
I have already my DM shell, however I have not looked yet how make it in the context of MG. Now is probably time to do that.

DM shell
http://userweb.eng.gla.ac.uk/lukasz.kaczmarczyk/MoFem/html/group__dm.html

Regards,
Lukasz

--
What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
-- Norbert Wiener

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160315/0a2a9b7b/attachment-0001.html>