[petsc-dev] EM solver in XGC
Mark Adams
mfadams at lbl.gov
Tue Apr 14 13:46:03 CDT 2015
I have tried
call
MatConvert(a_ts%FJacobian,MATMPIAIJ,MAT_INITIAL_MATRIX,a_ts%FJacobian2,ierr);CHKERRQ(ierr)
and
call
MatConvert(a_ts%FJacobian,MATAIJ,MAT_INITIAL_MATRIX,a_ts%FJacobian2,ierr);CHKERRQ(ierr)
in 'next' and it is not finding the specialized routine.
I will look into it.
Mark
On Tue, Apr 14, 2015 at 1:41 PM, Matthew Knepley <knepley at gmail.com> wrote:
> On Tue, Apr 14, 2015 at 12:40 PM, Matthew Knepley <knepley at gmail.com>
> wrote:
>
>> On Tue, Apr 14, 2015 at 12:33 PM, Mark Adams <mfadams at lbl.gov> wrote:
>>
>>> MatNest does not support convert because it does not have a getrow.
>>>
>>
>> I think this got fixed in dev. Anyhow, this supports my contention that
>> MatNest is a piece of shit.
>>
>
> Use 'next', it has
>
> PETSC_EXTERN PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType
> newtype,MatReuse reuse,Mat *newmat)
>
> so it should work.
>
> Matt
>
>
>>
>>
>>> Should we use a 4 block fieldsplit?
>>>
>>
>> We give up our ability to check using exact solves, but go ahead and try
>> it if you think it will make
>> the physicists happy for now.
>>
>> Thanks,
>>
>> Matt
>>
>>
>>> On Tue, Apr 14, 2015 at 1:26 PM, Mark Adams <mfadams at lbl.gov> wrote:
>>>
>>>> OK,
>>>>
>>>> On Tue, Apr 14, 2015 at 1:25 PM, Matthew Knepley <knepley at gmail.com>
>>>> wrote:
>>>>
>>>>> On Tue, Apr 14, 2015 at 12:23 PM, Mark Adams <mfadams at lbl.gov> wrote:
>>>>>
>>>>>> I moved this to after the setup, which it has to do because it does
>>>>>> not have the prefix or set options.
>>>>>>
>>>>>> It is not failing with not finding the IS. Will MatNest be able to
>>>>>> merge two matrices together like this?
>>>>>>
>>>>>> I create a 4x4 MatNext but have just two fields in fieldSplit.
>>>>>>
>>>>>
>>>>> Yep, I believe that will work. If not, just convert the MatNest to an
>>>>> AIJ for now with MatConvert().
>>>>>
>>>>> Matt
>>>>>
>>>>>
>>>>>> On Tue, Apr 14, 2015 at 1:18 PM, Matthew Knepley <knepley at gmail.com>
>>>>>> wrote:
>>>>>>
>>>>>>> On Tue, Apr 14, 2015 at 12:13 PM, Mark Adams <mfadams at lbl.gov>
>>>>>>> wrote:
>>>>>>>
>>>>>>>> I get an error where it looks like the number fields in field split
>>>>>>>> is < 2. Do I need to set the size?
>>>>>>>>
>>>>>>>
>>>>>>> No, each call to SetIS() just appends to the field list.
>>>>>>>
>>>>>>>
>>>>>>>> I use this code but I do this before I have setup the TS. Should I
>>>>>>>> do this after I setup the TS?
>>>>>>>>
>>>>>>>
>>>>>>> This is why I hate using the API since it is difficult to get all
>>>>>>> the setup done at the right time. First,
>>>>>>> when you call PCFieldSplitSetIS(), is the PC type set?
>>>>>>>
>>>>>>> Thanks,
>>>>>>>
>>>>>>> Matt
>>>>>>>
>>>>>>>
>>>>>>>> ! setup field split
>>>>>>>> call TSGetSNES(a_ts,snes,ierr);CHKERRQ(ierr)
>>>>>>>> call SNESGetKSP(snes,ksp,ierr);CHKERRQ(ierr)
>>>>>>>> call KSPGetPC(ksp,pc,ierr);CHKERRQ(ierr)
>>>>>>>>
>>>>>>>> call
>>>>>>>> MatGetSize(a_ts%mass,PETSC_NULL_INTEGER,npetsc,ierr);CHKERRQ(ierr)
>>>>>>>> npetsc = 2*npetsc ! two blocks in each partition
>>>>>>>> j = 0
>>>>>>>> call ISCreateStride(w_comm,npetsc,j,ione,is,ierr);CHKERRQ(ierr)
>>>>>>>> call PCFieldSplitSetIS(pc,'dt',is,ierr);CHKERRQ(ierr)
>>>>>>>> call ISDestroy(is,ierr);CHKERRQ(ierr)
>>>>>>>> j = npetsc
>>>>>>>> call ISCreateStride(w_comm,npetsc,j,ione,is,ierr);CHKERRQ(ierr)
>>>>>>>> call PCFieldSplitSetIS(pc,'aux',is,ierr);CHKERRQ(ierr)
>>>>>>>> call ISDestroy(is,ierr);CHKERRQ(ierr)
>>>>>>>>
>>>>>>>>
>>>>>>>> On Tue, Apr 14, 2015 at 12:24 PM, Matthew Knepley <
>>>>>>>> knepley at gmail.com> wrote:
>>>>>>>>
>>>>>>>>> On Tue, Apr 14, 2015 at 11:16 AM, Mark Adams <mfadams at lbl.gov>
>>>>>>>>> wrote:
>>>>>>>>>
>>>>>>>>>> Should I just use
>>>>>>>>>>
>>>>>>>>>> PetscErrorCode VecGetSubVector(Vec X,IS is,Vec *Y)
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>> Yep, that is exactly what I use it for.
>>>>>>>>>
>>>>>>>>> Matt
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>> ?
>>>>>>>>>>
>>>>>>>>>> On Tue, Apr 14, 2015 at 12:15 PM, Mark Adams <mfadams at lbl.gov>
>>>>>>>>>> wrote:
>>>>>>>>>>
>>>>>>>>>>> I am thinking that I should create an IS that is my local part
>>>>>>>>>>> of each field and then get a "sub vector" with that somehow. Looking into
>>>>>>>>>>> this now.
>>>>>>>>>>>
>>>>>>>>>>> On Tue, Apr 14, 2015 at 12:00 PM, Mark Adams <mfadams at lbl.gov>
>>>>>>>>>>> wrote:
>>>>>>>>>>>
>>>>>>>>>>>> Whoop, old code still there.
>>>>>>>>>>>>
>>>>>>>>>>>> I now have to get access to vectors somehow.
>>>>>>>>>>>>
>>>>>>>>>>>> I now get vectors for a field and do stuff with them like:
>>>>>>>>>>>>
>>>>>>>>>>>> subroutine FormIFunction(ts,t_dummy,a_XX,a_Xdot,a_FF,a_ts,ierr)
>>>>>>>>>>>>
>>>>>>>>>>>> call TSGetDM(ts,dm,ierr);CHKERRQ(ierr)
>>>>>>>>>>>>
>>>>>>>>>>>> call
>>>>>>>>>>>> DMCompositeGetAccessArray(dm,a_FF,ifour,PETSC_NULL_INTEGER,Fsub,ierr);CHKERRQ(ierr)
>>>>>>>>>>>> call
>>>>>>>>>>>> DMCompositeGetAccessArray(dm,a_Xdot,ifour,PETSC_NULL_INTEGER,Xdotsub,ierr);CHKERRQ(ierr)
>>>>>>>>>>>> one = 1d0
>>>>>>>>>>>> call VecAXPY(Fsub(1),one,Xdotsub(1),ierr);CHKERRQ(ierr)
>>>>>>>>>>>> call VecAXPY(Fsub(2),one,Xdotsub(2),ierr);CHKERRQ(ierr)
>>>>>>>>>>>> call
>>>>>>>>>>>> DMCompositeRestoreAccessArray(dm,a_Xdot,ifour,PETSC_NULL_INTEGER,Xdotsub,ierr);CHKERRQ(ierr)
>>>>>>>>>>>> call
>>>>>>>>>>>> DMCompositeRestoreAccessArray(dm,a_FF,ifour,PETSC_NULL_INTEGER,Fsub,ierr);CHKERRQ(ierr)
>>>>>>>>>>>>
>>>>>>>>>>>> call
>>>>>>>>>>>> MatMultAdd(a_ts%FJacobian,a_XX,a_FF,a_FF,ierr);CHKERRQ(ierr)
>>>>>>>>>>>>
>>>>>>>>>>>> I need to get rid of DMCompositeGetAccessArray. What do I do?
>>>>>>>>>>>>
>>>>>>>>>>>> Thanks,
>>>>>>>>>>>> Mark
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>> On Tue, Apr 14, 2015 at 11:36 AM, Matthew Knepley <
>>>>>>>>>>>> knepley at gmail.com> wrote:
>>>>>>>>>>>>
>>>>>>>>>>>>> On Tue, Apr 14, 2015 at 10:35 AM, Mark Adams <mfadams at lbl.gov>
>>>>>>>>>>>>> wrote:
>>>>>>>>>>>>>
>>>>>>>>>>>>>> And is this the solver (appended).
>>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>> Do not call SNESSetDM(). We are not using DM at all in this
>>>>>>>>>>>>> code.
>>>>>>>>>>>>>
>>>>>>>>>>>>> Matt
>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>>> I am getting an error at this last line:
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM()
>>>>>>>>>>>>>> @*/
>>>>>>>>>>>>>> PetscErrorCode SNESSetDM(SNES snes,DM dm)
>>>>>>>>>>>>>> {
>>>>>>>>>>>>>> PetscErrorCode ierr;
>>>>>>>>>>>>>> KSP ksp;
>>>>>>>>>>>>>> DMSNES sdm;
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> PetscFunctionBegin;
>>>>>>>>>>>>>> PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
>>>>>>>>>>>>>> if (dm) {ierr =
>>>>>>>>>>>>>> PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);}
>>>>>>>>>>>>>> if (snes->dm) { /* Move the DMSNES context
>>>>>>>>>>>>>> over to the new DM unless the new DM already has one */
>>>>>>>>>>>>>> if (snes->dm->dmsnes && snes->dmAuto && !dm->dmsnes) {
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> The problem is that the argument 'dm' to this method is NULL
>>>>>>>>>>>>>> so I get a segv at the last term of this conditional.
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> Mark
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> -adv_ts_type beuler
>>>>>>>>>>>>>> -adv_ts_monitor
>>>>>>>>>>>>>> -adv_ts_view
>>>>>>>>>>>>>> -adv_ts_max_step 1
>>>>>>>>>>>>>> -adv_snes_max_it 1 # make >> 1 for
>>>>>>>>>>>>>> nonlinear solve
>>>>>>>>>>>>>> -adv_snes_lag_preconditioner -1 # make > 0 to update PC
>>>>>>>>>>>>>> -adv_snes_lag_jacobian -1 # make > 0 for
>>>>>>>>>>>>>> nonlinear Jacobian
>>>>>>>>>>>>>> -adv_snes_monitor
>>>>>>>>>>>>>> -adv_snes_converged_reason
>>>>>>>>>>>>>> -adv_ksp_monitor
>>>>>>>>>>>>>> -adv_ksp_converged_reason
>>>>>>>>>>>>>> -adv_ksp_type richardson
>>>>>>>>>>>>>> -adv_pc_type fieldsplit
>>>>>>>>>>>>>> #-adv_pc_fieldsplit_type multiplicative
>>>>>>>>>>>>>> -adv_pc_fieldsplit_type schur
>>>>>>>>>>>>>> -adv_pc_fieldsplit_schur_fact_type full
>>>>>>>>>>>>>> -adv_pc_fieldsplit_schur_precondition full
>>>>>>>>>>>>>> #-adv_pc_fieldsplit_0_fields 0,1
>>>>>>>>>>>>>> #-adv_pc_fieldsplit_1_fields 2,3
>>>>>>>>>>>>>> -adv_fieldsplit_dt_pc_type lu
>>>>>>>>>>>>>> -adv_fieldsplit_aux_pc_type lu
>>>>>>>>>>>>>> -adv_fieldsplit_dt_pc_factor_mat_solver_package superlu
>>>>>>>>>>>>>> -adv_fieldsplit_aux_pc_factor_mat_solver_package superlu
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>> --
>>>>>>>>>>>>> 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
>>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> --
>>>>>>>>> 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
>>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> --
>>>>>>> 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
>>>>>>>
>>>>>>
>>>>>>
>>>>>
>>>>>
>>>>> --
>>>>> 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
>>>>>
>>>>
>>>>
>>>
>>
>>
>> --
>> 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
>>
>
>
>
> --
> 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-dev/attachments/20150414/df4a8a7c/attachment.html>
More information about the petsc-dev
mailing list