# Problem with vector fespace and result of matrix inversion

**URL:** <https://community.freefem.org/t/problem-with-vector-fespace-and-result-of-matrix-inversion/2296>\
**Category:** General Discussion\
**Created:** [February 27, 2023, 10:51pm UTC](https://community.freefem.org/t/problem-with-vector-fespace-and-result-of-matrix-inversion/2296 "2023-02-27T22:51:37Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [February 27, 2023, 10:51pm UTC](https://community.freefem.org/t/problem-with-vector-fespace-and-result-of-matrix-inversion/2296/1 "2023-02-27T22:51:37Z")

</div>

I copied one of the mpi heat examples and tried to adapt it into a  
multi-species diffusion problem based on the syntax from other examples.  
I can get the A^-1 \* b to execute but using the result is huge problem.  
In the past, I had just used scalar varf functions and made my own composite  
block diagonal matrix and it inverted and multipied just fine. I thought this  
would make better use of ff features but can’t figure out what is wrong.  
Now it complains during plaotting tht the solution field is empty except for the  
first vector comooent.  
Thanks.

```auto
[gluefemafk.edp|attachment](upload://6KxuL1CzyzrDkXuAcRHUxcirpOx.edp) (7.3 KB)

```

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [March 1, 2023, 12:25pm UTC](https://community.freefem.org/t/problem-with-vector-fespace-and-result-of-matrix-inversion/2296/2 "2023-03-01T12:25:35Z")

</div>

It looks like my code upload in the last message failed but I was able to  
get the code to run just by using “problem” although it is only using  
one thread. I realized in my home brew FD code it was worth making  
each species diffuse on its own core and that is only a few lines of  
openmp code. I guess the matrix is still block diagonal with coupling  
on the RHS as reactions occur with the old solution values.  
I’m not claiming the output is correct or useful just that it runs and gives  
plausible output. The mesh is designed to be compatible with  
my home brew FD code 🙂 If I do anything serious I will probably try to  
also run by FD code and maybe another FEM like libmesh for validation  
of my FF code. Thanks.

 ![femglue](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/e/e12b9890229197e5c4e321a8c21dadd859d80410.jpeg)

[gluefemvec.edp](https://community.freefem.org/uploads/short-url/4T0lz1CNiZPMbXYnoLnV47SZ5MW.edp) (5.2 KB)

```auto
 cat gluefemvec.edp // run with MPI: ff-mpirun -np 4 script.edp
// NBPROC 4
//
// code borrowed from examples/ffddm/heat-torus-3d-surf.edp

/*

@software{mjm_cpp_$PROJ,
  author = {Michael J Marchywka},
  title = {$PROJ},
abstract=(),
institution={},
license={Knowledge sir should be free to all },
publisher={Mike Marchywka},
email={marchywka@hotmail.com},
authorid={orcid.org/0000-0001-9237-455X},
  filename = {$PROJ},
  url = {},
  version = {0.0.0},
  date-started = {$DATE}
}

*/

// finally had to remove all that junk 

//load "hpddm" // HPDDM plugin
load "medit" // HPDDM plugin
//macro dimension()2// EOM // 2D, 3D, or 3S
include "macro_ddm.idp" // additional DDM functions

//macro def(i)i// EOM // scalar field definition

// syntax errors due to forgotten macro are a huge problem
// and not clear from eror message !!!!!!!!!!!!!!!!!!!!
//macro init(i)i// EOM // scalar field initialization
macro Grad(u)[dx(u), dy(u)]// EOM // two-dimensional gradient
func Pk = P2; // finite element space

int s = getARGV("-split", 1); // refinement factor
int sz = getARGV("-sz", 50); // refinement factor
int sz2 = getARGV("-sz2", 10); // refinement factor
real dt = getARGV("-dt", 1e-4); // time step
int iMax = getARGV("-iMax", 50); // number of iterations

int nx =getARGV("-nx", 30);
int ny = getARGV("-ny", 30);
cout<<" sz "<<sz<<" sz2 "<<sz2<<endl;
//meshS Th = square3(nx, ny, [(R+r*cos(2*pi*x))*cos(2*pi*y), (R+r*cos(2*pi*x))*sin(2*pi*y), r*sin(2*pi*x)]);
/////////////////////////////////////////////////
load "medit"
//https://doc.freefem.org/documentation/BEM.html
load "msh3"
// 2real R = 3, r=1, h=0.2;
func real a(real t) { real slow=.3*1.5; real fast=1.5*1.5; real bp1=.5-1.2*sz2/sz; real bp2=1.0-bp1;real norm=fast*(bp1+1-bp2)+slow*(bp2-bp1); if (t<bp1) return fast*t/norm; if ( t<bp2) return bp1*fast/norm+(t-bp1)*slow/norm; return bp1*fast/norm+(bp2-bp1)*slow/norm+(t-bp2)*fast/norm; } 

func fx= -sz + a(x)*2*sz;
func fy= -sz + a(y)*2*sz;
x=1;
cout<<" mnx mny "<<nx<<" "<<ny<<" "<<fx<<endl;
mesh Th = square(nx,ny,[fx,fy]);

//if (mpirank==0) 
medit("asdf",Th);

//////////////////////////////////////////////////
macro WVDIFF(u,v,w)(((u*Grad(w)'*Grad(v) ))) // EOM

fespace Wh(Th, [Pk,Pk,Pk,Pk,Pk,Pk]); // local finite element space
fespace Vh(Th, Pk); // local finite element space
int[int][int] intersection; // local-to-neighbors renumbering
real[int] Du; // partition of unity
//buildMinimalist(Th, intersection, Du, Pk)
//buildMinimalist(Th, intersection, Du, [Pk,Pk,Pk,Pk,Pk,Pk])

matrix<real> mat; // local operator
real f1z= getARGV("-f1z", 100); //
real f1az= getARGV("-f1az", 0); //
real f2z= getARGV("-f2z", 100); //
real f2az= getARGV("-f2az", 0); //
real fgnz= getARGV("-fgnz", 100); //
real fbz= getARGV("-fbz", 0); //
real tfz= getARGV("-tfz", 20); //
real c2= getARGV("-c2", .001); //
real cfgn= getARGV("-cfgn", .001); //

real df1= getARGV("-df1", 1); //
real df1a= getARGV("-df1a", .001); //
real df2= getARGV("-df2", 1); //
real df2a= getARGV("-df2a", .0001); //
real dfgn= getARGV("-dfgn", 1); //
real dfb= getARGV("-dfb", .001); //
real dzed= getARGV("-dzed",6000); //
real edfb= getARGV("-edfb",100); //
real cdfb= getARGV("-cdfb",.1); //

func real tfx( real xx, real yy)
{
//return 1000;
//cout<<xx<<" "<<sz2<<endl;
if ((xx<=sz2)&&(xx>= -sz2))  
if ((yy<=sz2)&&(yy>= -sz2)) return tfz;
return 0;
} 
Vh tf= tfx(x,y);
Wh<real> [f1,f1a,f2,f2a,fgn,fb]=[f1z,f1az,f2z,f2az,fgnz,fbz];
Wh<real> [f1w,f1aw,f2w,f2aw,fgnw,fbw]=[f1z,f1az,f2z,f2az,fgnz,fbz];
Wh<real> [vf1,vf1a,vf2,vf2a,vfgn,vfb];
Vh<real> D= dzed/(1+cdfb*exp(fb*edfb)) ;
problem vDiff([f1,f1a,f2,f2a,fgn,fb],
 [vf1,vf1a,vf2,vf2a,vfgn,vfb],solver=GMRES,init=1) =
//problem vDiff(f1,f1a,f2,f2a,fgn,fb, vf1,vf1a,vf2,vf2a,vfgn,vfb,solver=CG,init=j) =
int2d(Th) (
+dt*(df1*WVDIFF(D,f1,vf1)+df1a*WVDIFF(D,f1a,vf1a)+df2*WVDIFF(D,f2,vf2)  
+df2a*WVDIFF(D,f2a,vf2a)+dfgn*WVDIFF(D,fgn,vfgn)+dfb*WVDIFF(D,fb,vfb)  
)
+(f1*vf1+f1a*vf1a+f2*vf2+f2a*vf2a+fgn*vfgn+fb*vfb)
 )
-int2d(Th)(f1w*vf1+f1aw*vf1a+f2w*vf2+f2aw*vf2a+fgnw*vfgn+fbw*vfb)
+int2d(Th)(dt*f1w*vf1*tf)
+int2d(Th)(-dt*f1w*vf1a*tf)
+int2d(Th)(dt*f2w*vf2*f1aw*c2)
+int2d(Th)(-dt*f2w*vf2a*f1aw*c2)
+int2d(Th)(dt*fgnw*vfgn*f2aw*cfgn)
+int2d(Th)(-dt*fgnw*vfb*f2aw*cfgn)
//+on(1,2,3,4,[f1,f1a,f2,f2a,fgn,fb]=[f1z,f1az,f2z,f2az,fgnz,fbz])
+on(1,2,3,4,f1=f1z) +on(1,2,3,4,f1a=f1az) +on(1,2,3,4,f2=f2z)
+on(1,2,3,4,f2a=f2az) +on(1,2,3,4,fgn=fgnz) +on(1,2,3,4,fb=fbz)
;
for(int i=0; i<10000; ++i)
{
vDiff;
[f1w,f1aw,f2w,f2aw,fgnw,fbw]=[f1,f1a,f2,f2a,fgn,fb];
//medit("asdf",Th,fb);
  macro myplot()cmm = "Global solution at iteration " + i, fill = 1, value = 1//

plot(fb,myplot);
if ((i%1000)==9){
 medit("f1",Th,f1);
 medit("f1a",Th,f1a);
 medit("f2",Th,f2);
 medit("f2a",Th,f2a);
 medit("fgn",Th,fgn);
 medit("fb",Th,fb);
 medit("D",Th,D);
}

 D= dzed/(1+cdfb*exp(fb*edfb)) ;
} // i 

```