1 /* this sets up the vectors to call Fastcap's ComputePsi */
2 /* It will be called twice for each coordinate direction. Once for real
3 and once for imaginary */
4
5 #include "induct.h"
6
7 /* Vs will contain the result, Im is the 'q', Size is the size of vectors. */
8 /* This will alter Im. Im = Precond*Im */
SetupComputePsi(Vs,sys,Im,size,chglist,w,R,indsys)9 SetupComputePsi(Vs, sys, Im, size, chglist, w, R, indsys)
10 CX *Vs, *Im;
11 int size;
12 ssystem *sys;
13 charge *chglist;
14 double w; /* radian frequency */
15 double *R; /* resistance vector */
16 SYS *indsys;
17 {
18 extern double dirtime;
19 double *q, *p;
20 static CX *Ib = NULL, *Vb = NULL, *Vdirect = NULL, *ctemp;
21 int branches;
22 CX temp;
23 MELEMENT *mtemp;
24 charge *chg;
25 int i, j;
26 double rtemp;
27 MELEMENT **Mtrans, **Mlist;
28 double maxdiff,pdiff;
29 int maxindx;
30 int ind_opcnt_mult = 0, ind_opcnt_real = 0;
31
32 branches = indsys->num_fils;
33 Mtrans = indsys->Mtrans;
34 Mlist = indsys->Mlist;
35
36 if (Ib == NULL) {
37 Ib = (CX *)MattAlloc(branches, sizeof(CX));
38 Vb = (CX *)MattAlloc(branches, sizeof(CX));
39 ctemp = (CX *)MattAlloc(size, sizeof(CX));
40 #ifndef NODEBUG
41 Vdirect = (CX *)MattAlloc(branches, sizeof(CX));
42 #endif
43 }
44
45 for(i = 0; i < branches; i++)
46 Vb[i] = CXZERO;
47
48 q = sys->q;
49 p = sys->p;
50 ASSERT(size == indsys->num_mesh);
51
52 if (indsys->precond_type == LOC) {
53 multPrecond(indsys->Precond, Im, ctemp, size);
54 for(i = 0; i < size; i++)
55 Im[i] = ctemp[i];
56 }
57 else if (indsys->precond_type == SPARSE)
58 spSolve(indsys->sparMatrix, Im, Im);
59
60 /* do Ib = Mtrans*Im */
61 for(i = 0; i < branches; i++) {
62 Ib[i] = CXZERO;
63 for(mtemp = Mtrans[i]; mtemp != NULL; mtemp = mtemp->mnext) {
64 if (mtemp->sign == 1)
65 cx_add(Ib[i], Ib[i], Im[mtemp->filindex]);
66 else
67 cx_sub(Ib[i], Ib[i], Im[mtemp->filindex]);
68 }
69 }
70
71 /* Evaluate M*L*Mt*Im = M*L*Ib using the multipole algorithm */
72
73 /* Do all of the non-direct parts first */
74 sys->DirectEval = FALSE;
75 for(i = 0; i < 3; i++) { /* for each of the coordinate directions */
76
77 /* do the real part */
78 for(chg = chglist; chg != NULL; chg = chg->next) {
79 /* fill the pseudo-charge vector */
80 q[chg->index] = Ib[chg->fil->filnumber].real*chg->fil->lenvect[i];
81 #if OPCNT == ON
82 ind_opcnt_mult++;
83 #endif
84
85 }
86 computePsi(sys, q, p, branches, chglist);
87 for(chg = chglist; chg != NULL; chg = chg->next) {
88 /* add potential due to i direction */
89 Vb[chg->fil->filnumber].real += p[chg->index]*chg->fil->lenvect[i]*MUOVER4PI;
90 #if OPCNT == ON
91 ind_opcnt_mult++;
92 #endif
93 }
94
95 /* do the imaginary part */
96 for(chg = chglist; chg != NULL; chg = chg->next) {
97 /* fill the pseudo-charge vector */
98 q[chg->index] = Ib[chg->fil->filnumber].imag*chg->fil->lenvect[i];
99 #if OPCNT == ON
100 ind_opcnt_mult++;
101 #endif
102 }
103 computePsi(sys, q, p, branches, chglist);
104 for(chg = chglist; chg != NULL; chg = chg->next) {
105 /* add potential due to i direction */
106 Vb[chg->fil->filnumber].imag += p[chg->index]*chg->fil->lenvect[i]*MUOVER4PI;
107 #if OPCNT == ON
108 ind_opcnt_mult++;
109 #endif
110 }
111
112 }
113
114 /* do the direct parts */
115 sys->DirectEval = TRUE;
116
117 /* do the real part of the Direct part */
118 for(i = 1; i <= branches; i++)
119 p[i] = 0;
120 for(chg = chglist; chg != NULL; chg = chg->next)
121 /* fill the pseudo-charge vector */
122 q[chg->index] = Ib[chg->fil->filnumber].real;
123
124 /* starttimer; */
125 mulDirect(sys);
126 mulEval(sys);
127 /* stoptimer; */
128 dirtime += dtime;
129
130 for(chg = chglist; chg != NULL; chg = chg->next) {
131 /* add potential due to i direction */
132 Vb[chg->fil->filnumber].real += p[chg->index];
133 }
134
135 /* do the imaginary part of the Direct part */
136 for(i = 1; i <= branches; i++)
137 p[i] = 0;
138 for(chg = chglist; chg != NULL; chg = chg->next)
139 /* fill the pseudo-charge vector */
140 q[chg->index] = Ib[chg->fil->filnumber].imag;
141
142 /* starttimer; */
143 mulDirect(sys);
144 mulEval(sys);
145 /* stoptimer; */
146 dirtime += dtime;
147
148 for(chg = chglist; chg != NULL; chg = chg->next) {
149 /* add potential due to i direction */
150 Vb[chg->fil->filnumber].imag += p[chg->index];
151 }
152
153 /* do Vs = M*Vb*jw */
154 for(i = 0; i < size; i++) {
155 Vs[i] = CXZERO;
156 for(mtemp = Mlist[i]; mtemp != NULL; mtemp = mtemp->mnext)
157 if (mtemp->sign == 1)
158 cx_add(Vs[i], Vs[i], Vb[mtemp->filindex]);
159 else
160 cx_sub(Vs[i], Vs[i], Vb[mtemp->filindex]);
161
162 /* multiply by jw */
163 rtemp = -Vs[i].imag*w;
164 Vs[i].imag = Vs[i].real*w;
165 Vs[i].real = rtemp;
166 }
167
168 /* add in M*R*Mt*Im = M*R*Ib */
169 for(i = 0; i < size; i++) {
170 for(mtemp = Mlist[i]; mtemp != NULL; mtemp = mtemp->mnext) {
171 cx_scalar_mult(temp, mtemp->sign*R[mtemp->filindex], Ib[mtemp->filindex]);
172 cx_add(Vs[i], Vs[i], temp);
173 #if OPCNT == ON
174 ind_opcnt_mult+=2;
175 ind_opcnt_real+=2;
176 #endif
177 }
178 }
179
180 #if OPCNT == ON
181 printf("Inductance (mesh to branch) mults: %d\n",ind_opcnt_mult);
182 printf("Just doing MRMtIm: %d\n",ind_opcnt_real);
183 printops();
184 exit(0);
185 #endif
186
187 #ifdef NODEBUG
188 /* for debugging, compare to direct Vb = ZM Ib */
189 realmatCXvec(Vdirect, indsys->Z, Ib, branches);
190 maxdiff = 0;
191 maxindx = 0;
192 for(i = 0; i < branches; i++) {
193 if (cx_abs(Vb[i]) > 1e-23) {
194 cx_sub(temp, Vdirect[i], Vb[i]);
195 pdiff = cx_abs(temp )/cx_abs(Vb[i]) ;
196 }
197 else
198 pdiff = cx_abs(Vb[i]);
199
200 if (pdiff > maxdiff) {
201 maxdiff = pdiff;
202 maxindx = i;
203 }
204 }
205 if (maxdiff < .3)
206 printf("maxdiff: %g Vb[%d]=%g Vdirect[%d]=%g\n",
207 maxdiff,maxindx,cx_abs(Vb[maxindx]),maxindx,cx_abs(Vdirect[maxindx]));
208 else
209 printf("***maxdiff: %g Vb[%d]=%g Vdirect[%d]=%g***\n",
210 maxdiff,maxindx,cx_abs(Vb[maxindx]),maxindx,cx_abs(Vdirect[maxindx]));
211
212
213 #endif
214 }
215
realmatCXvec(y,A,x,size)216 realmatCXvec(y, A, x, size)
217 CX *y, *x;
218 double **A;
219 int size;
220 {
221 int i, j;
222 CX temp;
223
224 for (i = 0; i < size; i++) {
225 y[i] = CXZERO;
226 for(j = 0; j < size; j++) {
227 cx_scalar_mult(temp, A[i][j], x[j]);
228 cx_add(y[i], y[i], temp);
229 }
230 }
231 }
232
233 /* this function fixes Eval matrices which are computed directly */
234 /* This is necessary since direct mutual terms are not componentwise,
235 but the multipole routines are called once for each component direction.
236 Basically, componentwise multiplication will cause the elements
237 to be multiplied by the dot product of the fil->lenvect vectors of
238 the two filaments. This will divide that product out. Also, MUOVER4PI
239 must also be divided out
240 */
241
fixEvalDirect(qchgs,numqchgs,is_dummy,pchgs,numpchgs,mat)242 fixEvalDirect(qchgs, numqchgs, is_dummy, pchgs, numpchgs, mat)
243 charge **qchgs, **pchgs;
244 int numqchgs, numpchgs;
245 int *is_dummy;
246 double **mat;
247 {
248 int i,j, k;
249 double dotprod, magi, magj;
250 double *lenvecti, *lenvectj;
251 static double eps = EPS;
252
253 for(i = 0; i < numpchgs; i++) {
254 lenvecti = pchgs[i]->fil->lenvect;
255 magi = 0;
256 for(k = 0; k < 3; k++)
257 magi += lenvecti[k]*lenvecti[k];
258 for(j = 0; j < numqchgs; j++) {
259 lenvectj = qchgs[j]->fil->lenvect;
260 magj = dotprod = 0;
261 for(k = 0; k < 3; k++) {
262 magj += lenvectj[k]*lenvectj[k];
263 dotprod += lenvecti[k]*lenvectj[k];
264 }
265 if (fabs(dotprod)/sqrt(magi*magj) > EPS) /* filaments aren't perpendicular */
266 mat[i][j] = mat[i][j]/(dotprod*MUOVER4PI);
267 else { /* if they are, mat[i][j] == 0.0, hopefully */
268 if (mat[i][j] != 0.0)
269 printf("Warning: dot product = %lg < EPS, but mat[i][j] = %lg\n",dotprod, mat[i][j]);
270 }
271 }
272 }
273 }
274