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.
gluefemvec.edp (5.2 KB)
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
