<div dir="ltr">Thanks for your quick answer. <div><br></div><div>Yeah, the Vec I passed to KSPSolve is in the wrong blocksize. Only the rhs Vec is blocked, the solution Vec is not blocked.</div><div><br></div><div>Add  VecSetBlockSize before KSPSolve solved the problem. Thanks you very much.</div>
</div><div class="gmail_extra"><br><br><div class="gmail_quote">On Tue, Apr 29, 2014 at 4:06 PM, Jed Brown <span dir="ltr"><<a href="mailto:jed@jedbrown.org" target="_blank">jed@jedbrown.org</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div class="">Song Gao <<a href="mailto:song.gao2@mail.mcgill.ca">song.gao2@mail.mcgill.ca</a>> writes:<br>
<br>
> Dear Petsc users,<br>
><br>
> We are working on a matrix-free solver for NS equations. We first developed<br>
> the solver in the SNES framework. But because we want to take control of<br>
> the Newton iterations and reuse some of our existing code, we stop using<br>
> SNES and begin developing the matrix-free solver in KSP framework.<br>
><br>
> But we have an issue when using MFFD and VecSetValuesBlocked together.<br>
> My code looks like<br>
><br>
> call MatMFFDSetFunction(pet_mat_mf, flowsol_ng_mf,  ctx, ierpetsc)<br>
> call MatMFFDSetBase(pet_mat_mf, pet_solu_snes, PETSC_NULL_OBJECT, ierpetsc)<br>
><br>
> The vec pet_solu_snes are set to blocksize 5.<br>
<br>
</div>Did you set the block size for the vector passed to KSP?  MatMFFD does<br>
not create its own vector (and the base vector would be the wrong size<br>
for a non-square matrix anyway).<br>
<div class="HOEnZb"><div class="h5"><br>
> flowsol_ng_mf looks like:<br>
> subroutine flowsol_ng_mf ( ctxx, pet_solu_snesa, pet_rhs_snesa,  ierrpet )<br>
> ...............<br>
>     call VecGetBlockSize(pet_rhs_snesa,i,ierpetsc)<br>
>     write(*,*) "pet_rhs_snesa",i<br>
> .............<br>
>     call VecSetValuesBlocked ( pet_rhs_snesa, ndperl, pirow2, rhsloc,<br>
> ADD_VALUES, ierpetsc )<br>
> ..........<br>
> end<br>
><br>
> I'm worried about the block size of pet_rhs_snesa is not properly<br>
> initialized. So I checked the blocksize in the function evaluation.<br>
><br>
> The output with single CPU looks like:<br>
><br>
>   0 KSP preconditioned resid norm 7.818306841541e+00 true resid norm<br>
> 9.619278462343e-03 ||r(i)||/||b|| 1.000000000000e+00<br>
>  pet_rhs_snesa           1<br>
>  pet_rhs_snesa           1<br>
>  pet_rhs_snesa           5<br>
>   1 KSP preconditioned resid norm 1.739193723080e+00 true resid norm<br>
> 2.564460053330e-03 ||r(i)||/||b|| 2.665958848545e-01<br>
>  pet_rhs_snesa           1<br>
>  pet_rhs_snesa           5<br>
>   2 KSP preconditioned resid norm 3.816590845074e-01 true resid norm<br>
> 8.222646506071e-04 ||r(i)||/||b|| 8.548090730777e-02<br>
>  pet_rhs_snesa           1<br>
>  pet_rhs_snesa           5<br>
><br>
> I think the block size is not initialized correctly for the first call,<br>
> which leads to the misuse of  VecSetValuesBlocked.<br>
><br>
> One way to circumvent is to create a temporary vec with correct blocksize,<br>
> assemble it with VecSetValuesBlocked and then copy back. But Is there other<br>
> solutions? Thank you very much.<br>
><br>
><br>
><br>
> Song<br>
><br>
><br>
><br>
> On Sat, Mar 8, 2014 at 8:51 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>> wrote:<br>
><br>
>> On Sat, Mar 8, 2014 at 7:45 PM, Song Gao <<a href="mailto:song.gao2@mail.mcgill.ca">song.gao2@mail.mcgill.ca</a>> wrote:<br>
>><br>
>>> Thank you all.<br>
>>><br>
>>> Most of the difficulties come from the current structure of the program.<br>
>>> Its size and "age" make any radical modification a challenge. The large<br>
>>> uses of global variables and deprecated "goto" statements call for a<br>
>>> revision, which however is unlikely to occur, against our better<br>
>>> judgement...<br>
>>><br>
>>> That being said, the current tools provided byPETSc are sufficient,<br>
>>> however not necessarily convenient for our purposes. As a wishful thinking<br>
>>> we can say that the implementation of PETSc SNES features into existing<br>
>>> codes would be easier if the outer Newton-like loop could be managed by the<br>
>>> original code, out of the SNES context. On the other hand, we realize that<br>
>>> this might be in contrast with the requirements of PETSc itself.<br>
>>><br>
>><br>
>> I guess I wanted to know specifically why you cannot just dump all of the<br>
>> other things in the legacy code into SNESSetUpdate()<br>
>> or SNESSetConvergenceTest(). Was there something about the way we<br>
>> structured the hooks that made this impossible?<br>
>><br>
>>   Thanks,<br>
>><br>
>>      Matt<br>
>><br>
>><br>
>>> The application of the Eistenstat-Walker method together with a<br>
>>> Matrix-Free approach are of key importance. Thus, as kindly suggested by<br>
>>> Jed Brown, implementing our own version of EW scheme might turn out to be<br>
>>> the only way. We would be happy to provide more details if you think that<br>
>>> they might be helpful for the future developments of PETSc.<br>
>>><br>
>>><br>
>>> On Fri, Mar 7, 2014 at 2:39 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>>wrote:<br>
>>><br>
>>>> On Fri, Mar 7, 2014 at 1:21 PM, Song Gao <<a href="mailto:song.gao2@mail.mcgill.ca">song.gao2@mail.mcgill.ca</a>>wrote:<br>
>>>><br>
>>>>> Hello,<br>
>>>>><br>
>>>>> We are working on a legacy codes which solves the NS equations. The<br>
>>>>> codes take control of newton iterations itself and use KSP to solve the<br>
>>>>> linear system.<br>
>>>>><br>
>>>>> We modified the codes, used SNES to control the newton iteration and<br>
>>>>> changed the code to matrix free fashion. But the legacy codes did a lot of<br>
>>>>> other things between two newton iterations (such as output solution, update<br>
>>>>> variables....). I know we could use linesearchpostcheck but it is<br>
>>>>> difficulty to do that correctly. Therefore, we decide to go back to the KSP<br>
>>>>> framework but still use matrix free.<br>
>>>>><br>
>>>><br>
>>>> What makes it difficult?<br>
>>>><br>
>>>>   Thanks,<br>
>>>><br>
>>>>      Matt<br>
>>>><br>
>>>><br>
>>>>> When using SNES, we always use the runtime option -snes_ksp_ew, we<br>
>>>>> observe that for some test cases, the residual stalls without -snes_ksp_ew,<br>
>>>>> but converges with-snes_ksp_ew. So I'm thinking if it is possible to<br>
>>>>> use -snes_ksp_ew in KSP? Thanks in advance.<br>
>>>>><br>
>>>>><br>
>>>>><br>
>>>><br>
>>>><br>
>>>> --<br>
>>>> What most experimenters take for granted before they begin their<br>
>>>> experiments is infinitely more interesting than any results to which their<br>
>>>> experiments lead.<br>
>>>> -- Norbert Wiener<br>
>>>><br>
>>><br>
>>><br>
>><br>
>><br>
>> --<br>
>> What most experimenters take for granted before they begin their<br>
>> experiments is infinitely more interesting than any results to which their<br>
>> experiments lead.<br>
>> -- Norbert Wiener<br>
>><br>
</div></div></blockquote></div><br></div>