# Call PETSc matrix operator MatDiagonalSet (…) **URL:** \ **Category:** General Discussion\ **Created:** [June 21, 2021, 12:41pm UTC](https://community.freefem.org/t/call-petsc-matrix-operator-matdiagonalset/1050 "2021-06-21T12:41:33Z")\ **Posts on this page:** 19\ **Page:** 2
**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\ **Post date:** [June 28, 2021, 8:19pm UTC](https://community.freefem.org/t/call-petsc-matrix-operator-matdiagonalset/1050/25 "2021-06-28T20:19:50Z")
I’m not sure why that is. The parallel interpolation is not really designed for your problem (with meshes of greatly different dimensions). In your case, your leaflet is very small, so it may be counterproductive to distribute it across processes. Instead, you could, for example, let process #0 hold the solid mesh. Then, distribute the fluid mesh by making sure that the solid mesh only lives in the portion of the fluid mesh held by process #0. That way, you can do all fluid \<-\> solid interpolation on process #0, “in sequential” (which may be better since the domain is tiny), and do all fluid work in parallel. ---
**Author:** ![yongxing](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/yongxing/32/3632_2.png) [@yongxing](https://community.freefem.org/u/yongxing)\ **Post date:** [June 28, 2021, 8:48pm UTC](https://community.freefem.org/t/call-petsc-matrix-operator-matdiagonalset/1050/27 "2021-06-28T20:48:34Z")
Thank you, I will give it a try. However, this is just a test. The motivation is to simulate the whole heart including all the valves, in which case the solid (valves and tissues) would be large ---
**Author:** ![yongxing](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/yongxing/32/3632_2.png) [@yongxing](https://community.freefem.org/u/yongxing)\ **Post date:** [June 29, 2021, 4:00pm UTC](https://community.freefem.org/t/call-petsc-matrix-operator-matdiagonalset/1050/28 "2021-06-29T16:00:41Z")
It seems I always have to do the following code **together** (inside the time loop) in order to compute the interpolation matrix P ( **mesh Ths is moving but Th is static** ) createMat(Th, A, PV1)| createMat(Ths, B, PV1)| transferMat(Th, PV1, A, Ths, PV1, B, P) A = fluid(Rh, Rh); B = solid(Rhs,Rhs) One problem is I quickly ran out of memory because of (I guess) calling createMat() again and again. Another problem is matrix A is static, It’s better to put it outside the time loop. However, when I do this, it shows the following error: ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/1X/f30570a83503a472a57c948e8a78cf157e5cc76f.png) Can I somehow only **put the following two lines inside the time loop, and all the other three lines outside the time loop** because only Ths and matrix B change as the time involves ? transferMat(Th, PV1, A, Ths, PV1, B, P) B = solid(Rhs,Rhs) ---
**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\ **Post date:** [June 29, 2021, 7:24pm UTC](https://community.freefem.org/t/call-petsc-matrix-operator-matdiagonalset/1050/29 "2021-06-29T19:24:18Z")
Yes, of course you can. Maybe for good measure add `MatDestroy(P);` before the call to `transferMat`. ---
**Author:** ![yongxing](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/yongxing/32/3632_2.png) [@yongxing](https://community.freefem.org/u/yongxing)\ **Post date:** [July 1, 2021, 3:38pm UTC](https://community.freefem.org/t/call-petsc-matrix-operator-matdiagonalset/1050/31 "2021-07-01T15:38:08Z")
Dear Prof. Pierre Jolivet, Your know, after solving a PDE for **u** , we can use **dx(u)** to get the derivative of **u**. However, In the case parallel, it seems we cannot directly do so, because I found the results were different between using one processor and two processors. Do we have to gather **(centralise) u to one processor**, compute the derivative, and then distribute the results across different processors? or there is a simple trick? Best, Yongxing. ---
**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\ **Post date:** [July 1, 2021, 3:41pm UTC](https://community.freefem.org/t/call-petsc-matrix-operator-matdiagonalset/1050/32 "2021-07-01T15:41:47Z")
You can, but you need to synchronize the solution at ghost elements, see line 36 of [http://jolivet.perso.enseeiht.fr/FreeFem-tutorial/section\_8/example2.edp.html](http://jolivet.perso.enseeiht.fr/FreeFem-tutorial/section_8/example2.edp.html) (use dx instead of b). ---
**Author:** ![yongxing](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/yongxing/32/3632_2.png) [@yongxing](https://community.freefem.org/u/yongxing)\ **Post date:** [July 1, 2021, 7:03pm UTC](https://community.freefem.org/t/call-petsc-matrix-operator-matdiagonalset/1050/34 "2021-07-01T19:03:09Z")
Sorry, I am too silly to understand the details. I actually want to compute the deformation tensor **F** as follows: ```auto dispx[] += wx[]*dt; dispy[] += wy[]*dt; f11=dx(dispx)+1; f12=dy(dispx); f21=dx(dispy); f22=dy(dispy)+1; ``` However I found **F** was not right near the ghost elements. Do you mean I should put **exchange(A, dx, scaled = true); exchange(A, dy, scaled = true);** Between the above two lines of code? How should I choose matrix A? I have tried…which was problematic. ---
**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\ **Post date:** [July 1, 2021, 7:16pm UTC](https://community.freefem.org/t/call-petsc-matrix-operator-matdiagonalset/1050/35 "2021-07-01T19:16:40Z")
You need to `exchange` your f11, f12, etc. `A` must be a `Mat` defined using the same `fespace` as your f11, f12, etc. ---
**Author:** ![yongxing](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/yongxing/32/3632_2.png) [@yongxing](https://community.freefem.org/u/yongxing)\ **Post date:** [July 1, 2021, 8:04pm UTC](https://community.freefem.org/t/call-petsc-matrix-operator-matdiagonalset/1050/36 "2021-07-01T20:04:09Z")
Thank you vert much! I did the following: ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/1X/a83ac78427ed598f09394827a6ac7e4831361e3b.png) but it shows the following **operator error:** ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/1X/398097b3cc8908b55a1a62278a365f822d462846.png) ---
**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\ **Post date:** [July 1, 2021, 8:07pm UTC](https://community.freefem.org/t/call-petsc-matrix-operator-matdiagonalset/1050/37 "2021-07-01T20:07:31Z")
f11[] instead of f11… ---
**Author:** ![yongxing](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/yongxing/32/3632_2.png) [@yongxing](https://community.freefem.org/u/yongxing)\ **Post date:** [July 5, 2021, 9:19am UTC](https://community.freefem.org/t/call-petsc-matrix-operator-matdiagonalset/1050/39 "2021-07-05T09:19:55Z")
Dear prj, Can I ask you one more question about interpolation matrix: when I use transferMat(), the interpolation matrix it produces is always linear, or it uses the shape function, which could then be high order if using a high-order fespace? Best, Yongxing. ---
**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\ **Post date:** [July 5, 2021, 9:58am UTC](https://community.freefem.org/t/call-petsc-matrix-operator-matdiagonalset/1050/40 "2021-07-05T09:58:32Z")
It will use whatever you put in `PV1` and `P` (using your notations from [this](https://community.freefem.org/t/call-petsc-matrix-operator-matdiagonalset/1050/28) post), so it can be of high order. ---
**Author:** ![yongxing](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/yongxing/32/3632_2.png) [@yongxing](https://community.freefem.org/u/yongxing)\ **Post date:** [August 6, 2021, 9:30pm UTC](https://community.freefem.org/t/call-petsc-matrix-operator-matdiagonalset/1050/43 "2021-08-06T21:30:05Z")
Should I do **exchange (A, f11[], scaled=true)** immediately after **f11=dx(dispx)** or I can do operations on f11, such as computing the inverse of **F^{-1}=P** , even do some matrix operations such as **F\*F’+I=P** … then finally call **exchange (A, p11[], scaled=true)** or the order does not matter? ---
**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\ **Post date:** [August 7, 2021, 5:48am UTC](https://community.freefem.org/t/call-petsc-matrix-operator-matdiagonalset/1050/44 "2021-08-07T05:48:30Z")
How would that matter if none of the operations involve neither `f11` nor `A`? ---
**Author:** ![yongxing](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/yongxing/32/3632_2.png) [@yongxing](https://community.freefem.org/u/yongxing)\ **Post date:** [August 7, 2021, 7:04pm UTC](https://community.freefem.org/t/call-petsc-matrix-operator-matdiagonalset/1050/45 "2021-08-07T19:04:08Z")
Sorry, the operation does involve f11, which is the component of matrix F. In this case, should I perform the operation first, or exchange() first? ---
**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\ **Post date:** [August 7, 2021, 7:43pm UTC](https://community.freefem.org/t/call-petsc-matrix-operator-matdiagonalset/1050/46 "2021-08-07T19:43:25Z")
Well, if you want the proper value, of course you need to `exchange()` beforehand. I’m sorry, I’m not sure I understand the question correctly… ---
**Author:** ![yongxing](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/yongxing/32/3632_2.png) [@yongxing](https://community.freefem.org/u/yongxing)\ **Post date:** [August 7, 2021, 8:22pm UTC](https://community.freefem.org/t/call-petsc-matrix-operator-matdiagonalset/1050/47 "2021-08-07T20:22:54Z")
Some very simple examples are: a=dx(disp) b=a+1 c=a\*a d=sin(a) should I exchange ( … a…) first, then compute b, c or d? can I first compute b, c or d, then exchange (…b…), exchange (…c…) or exchange (…d…)? or they are the same? ---
**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\ **Post date:** [August 7, 2021, 8:34pm UTC](https://community.freefem.org/t/call-petsc-matrix-operator-matdiagonalset/1050/48 "2021-08-07T20:34:04Z")
I think it’s important you try to understand why it would be the same for `b` and `d` but different for `c` if you don’t do the `exchange()` beforehand. ---
**Author:** ![WeiQLiu](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/weiqliu/32/1043_2.png) [@WeiQLiu](https://community.freefem.org/u/WeiQLiu)\ **Post date:** [June 3, 2024, 3:09pm UTC](https://community.freefem.org/t/call-petsc-matrix-operator-matdiagonalset/1050/49 "2024-06-03T15:09:12Z")
Dear prj, I think your method is very interesting. That is, let process #0 hold the solid mesh. Then, distribute the fluid mesh by making sure that the solid mesh only lives in the portion of the fluid mesh held by process #0. While my specific physics problem may differ from those discussed here, I ecounter a similar issue regarding the coupling between small computational domain and large computational domain, for example, linear stability analysis of a single bubble moving in uniform flow. However, I don’t understand how exactly this method is implemented, for example, how to specify a specific domain on the partition mesh of process #0 in large computational domain？ Best, wqLiu [Previous page](https://community.freefem.org/t/call-petsc-matrix-operator-matdiagonalset/1050.md?page=1)