root/trunk/rheoplast/chternary.c

Revision 1461, 78.4 kB (checked in by hazelsct, 2 years ago)

Major overhaul, too extensive to comment. Includes start of Jacobian.
  • Property svn:keywords set to Author Date Id Revision
Line 
1/***************************************
2  $Header$
3
4  This provides a ternary Cahn-Hilliard module for Rheoplast.  It uses free
5  energy given by
6  +latex+\begin{equation}
7  +latex+  \label{cht_energy}
8  +latex+  {\cal F} = \int \left(G (C_2,C_3) +
9  +latex+    \frac{1}{2} \sum_{i,j=2,3}K_{ij}\nabla C_i \cdot \nabla C_j
10  +latex+  \right) dV,
11  +latex+\end{equation}
12  +latex+where $G(C_2,C_3)$ is the homogeneous free energy, whose derivatives
13  +latex+are described in sections \ref{func_psiprime2_chternary.c} and
14  +latex+\ref{func_psiprime3_chternary.c} (pages
15  +latex+\pageref{func_psiprime2_chternary.c} and
16  +latex+\pageref{func_psiprime3_chternary.c}).
17  +html+ the sum of homogeneous free energy (whose derivatives are described in
18  +html+ <a href="#func-psiprime2">psiprime2</a> and
19  +html+ <a href="#func-psiprime3">psiprime3</a>) and gradient penalty terms.
20  ***************************************/
21 
22
23#include "rheoplast.h"
24
25
26#undef __FUNCT__
27#define __FUNCT__ "pseudolog"
28
29/*++++++++++++++++++++++++++++++++++++++
30  Returns a "pseudo-log" value which doesn't FPE for values below 0.
31
32  inline PetscScalar pseudolog Return value.
33
34  PetscScalar arg Argument.
35  ++++++++++++++++++++++++++++++++++++++*/
36
37static inline PetscScalar pseudolog (PetscScalar arg)
38{ return (arg>1e-10) ? log(arg) : log(1e-10) + 1e10*(arg - 1e-10); }
39
40
41#undef __FUNCT__
42#define __FUNCT__ "psiprime2"
43
44/*++++++++++++++++++++++++++++++++++++++
45  This abstracts out the function for the derivative of homogeneous free energy
46  with respect to
47  +latex+$C_2$: $\partial G/\partial C_2 (C_2,C_3)$
48  +html+ <i>C</i><sub>2</sub>: <i>dG</i>/<i>dC</i><sub>2</sub>
49  +html+ (<i>C</i><sub>2</sub>, <i>C</i><sub>3</sub>)
50  (which is a component of the chemical potential of species 2) so it can be
51  easily modified.
52 
53  PetscScalar psiprime2 It returns the derivative of homogeneous free
54  energy.
55
56  PetscScalar C2 Concentration of species 2.
57
58  PetscScalar C3 Concentration of species 3.
59
60  chtparm *thechternary chternary parameter structure.
61  ++++++++++++++++++++++++++++++++++++++*/
62
63static inline PetscScalar psiprime2 (PetscScalar C2, PetscScalar C3,
64                                     chtparm *thechternary)
65{
66  int ierr;
67
68  if (thechternary->polymer)
69    {
70      /*+ In the polymer solution case,
71        +latex+$C_2$ and $C_3$
72        +html+ <i>C</i><sub>2</sub> and <i>C</i><sub>3</sub>
73        are solvent and polymer volume fractions
74        +latex+$\phi_s$ and $\phi_p$, and
75        +html+ <i>phi<sub>s</sub></i> and <i>phi<sub>p</sub></i>, and
76        the homogeneous free energy derivative is given by the Flory-Huggins
77        +latex+model:
78        +latex+\begin{eqnarray}
79        +latex+  \label{eq:psiprime2}
80        +latex+  \frac{\partial G}{\partial\phi_s} = & RT & \left(
81        +latex+    \frac{-1-\ln(1-\phi_s-\phi_p)}{V_n} +
82        +latex+    \frac{1+\ln \phi_s}{V_s} +
83        +latex+    \frac{\chi_{ns}(1-2\phi_s-\phi_p) +
84        +latex+      \frac{\partial\chi_{ns}}{\partial\phi_s}
85        +latex+      (1-\phi_s-\phi_p)\phi_s}{V_n} +
86        +latex+  \right. \\ \nonumber && \left.
87        +latex+    \frac{\chi_{sp}\phi_p +
88        +latex+      \frac{\partial\chi_{sp}}{\partial\phi_s}\phi_s\phi_p}{V_s}
89        +latex+    + \frac{-\chi_{np}\phi_p + \frac{\partial\chi_{np}}
90        +latex+      {\partial\phi_s}(1-\phi_s-\phi_p)\phi_p}{V_n}
91        +latex+  \right).
92        +latex+\end{eqnarray}
93        +html+ model.
94        +*/
95
96      /* PVDF-DMF-water parameters */
97      PetscScalar chi_ns = 0.218 + 0.276/(1-0.622*C2/(1-C3));
98      PetscScalar chi_sp = -1+0.5*C3;
99      PetscScalar chi_np = 3.5;
100
101      PetscScalar dchins_dC2 = 0.276/(1-0.622*C2/(1-C3))/(1-0.622*C2/(1-C3)) *
102        0.622*C2/(1-C3);
103      PetscScalar dchisp_dC2 = 0.0;
104      PetscScalar dchinp_dC2 = 0.0;
105
106      return ((-1.0-pseudolog(1.-C2-C3))/thechternary->V_n +
107              (1.+pseudolog(C2))/thechternary->V_s +
108              (chi_ns*(1.-2.*C2-C3) +
109               dchins_dC2*(1.-C2-C3)*C2)/thechternary->V_n +
110              (chi_sp*C3 + dchisp_dC2*C2*C3)/thechternary->V_s +
111              (-chi_np*C3 + dchinp_dC2*(1.-C2-C3)*C3)/thechternary->V_n) *
112        thechternary->RT;
113    }
114
115  /* Metal electrolysis parameters and homogeneous free energy */
116  PetscScalar chiMM =3.3;
117  PetscScalar chiMC =2.5;
118  PetscScalar chipMM =0.0;
119  PetscScalar chipMC =0.0;
120
121  return log(C2) - log(1.0-C2) + chiMM*C3*(1.0-2.*C2)-(C3-0.5) + chipMM*C3*(2.*C2-6.*C2*C2 +4.*C2*C2*C2);
122}
123
124
125#undef __FUNCT__
126#define __FUNCT__ "psiprime3"
127
128/*++++++++++++++++++++++++++++++++++++++
129  This abstracts out the function for the derivative of homogeneous free energy
130  with respect to
131  +latex+$C_3$: $\partial G/\partial C_3 (C_2,C_3)$
132  +html+ <i>C</i><sub>3</sub>: <i>dG</i>/<i>dC</i><sub>3</sub>
133  +html+ (<i>C</i><sub>2</sub>, <i>C</i><sub>3</sub>)
134  (which is a component of the chemical potential of species 3) so it can be
135  easily modified.
136 
137  PetscScalar psiprime3 It returns the derivative of homogeneous free
138  energy.
139
140  PetscScalar C2 Concentration of species 2.
141
142  PetscScalar C3 Concentration of species 3.
143
144  chtparm *thechternary chternary parameter structure.
145  ++++++++++++++++++++++++++++++++++++++*/
146
147static inline PetscScalar psiprime3 (PetscScalar C2, PetscScalar C3,
148                                     chtparm *thechternary)
149{
150  int ierr;
151
152  if (thechternary->polymer)
153    {
154      /*+ In the polymer solution case,
155        +latex+$C_2$ and $C_3$
156        +html+ <i>C</i><sub>2</sub> and <i>C</i><sub>3</sub>
157        are solvent and polymer volume fractions
158        +latex+$\phi_s$ and $\phi_p$, and
159        +html+ <i>phi<sub>s</sub></i> and <i>phi<sub>p</sub></i>, and
160        the homogeneous free energy derivative is given by the Flory-Huggins
161        +latex+model:
162        +latex+\begin{eqnarray}
163        +latex+  \label{eq:psiprime3}
164        +latex+  \frac{\partial G}{\partial\phi_p} = & RT & \left(
165        +latex+    \frac{-1-\ln(1-\phi_s-\phi_p)}{V_n} +
166        +latex+    \frac{1+\ln \phi_p}{V_p} +
167        +latex+    \frac{-\chi_{ns}\phi_s +
168        +latex+      \frac{\partial\chi_{ns}}{\partial\phi_p}
169        +latex+      (1-\phi_s-\phi_p)\phi_s}{V_n} +
170        +latex+  \right. \\ \nonumber && \left.
171        +latex+    \frac{\chi_{sp}\phi_s +
172        +latex+      \frac{\partial\chi_{sp}}{\partial\phi_p}\phi_s\phi_p}{V_s}
173        +latex+    + \frac{\chi_{np}(1-\phi_s-2\phi_p)+\frac{\partial\chi_{np}}
174        +latex+      {\partial\phi_p}(1-\phi_s-\phi_p)\phi_p}{V_n}
175        +latex+  \right).
176        +latex+\end{eqnarray}
177        +html+ model.
178        +*/
179
180      /* PVDF-DMF-water parameters */
181      PetscScalar chi_ns = 0.218 + 0.276/(1-0.622*C2/(1-C3));
182      PetscScalar chi_sp = -1+0.5*C3;
183      PetscScalar chi_np = 3.5;
184
185      PetscScalar dchins_dC3 = 0.276/ (1-0.622*C2/(1-C3)) / (1-0.622*C2/(1-C3))
186        * 0.622*C2/(1.-C3)/(1.-C3);
187      PetscScalar dchisp_dC3 = 0.5;
188      PetscScalar dchinp_dC3 = 0.0;
189
190      return ((-1.0-pseudolog(1.-C2-C3))/thechternary->V_n +
191              (1.+pseudolog(C3))/thechternary->V_p +
192              (-chi_ns*C2 + dchins_dC3*(1.-C2-C3)*C2)/thechternary->V_n +
193              ( chi_sp*C2 + dchisp_dC3*C2*C3)/thechternary->V_s +
194              ( chi_np*(1.-C2-2.*C3) + dchinp_dC3*(1.-C2-C3)*C3)
195              /thechternary->V_n) *
196        thechternary->RT;
197    }
198
199  PetscScalar chiMM =3.3;
200  PetscScalar chiMC =2.5;
201  PetscScalar chipMM =0.0;
202  PetscScalar chipMC =0.0;
203
204  return  chiMM*C2*(1.0-C2) - (C2-0.5) +log(C3)-log(1.0-C3) + chiMC*(1.0-2.*C3)+
205    chipMC*(2.*C3-6.*C3*C3 +4.*C3*C3*C3) + chipMM*(C2*C2)*(1-C2)*(1-C2);
206}
207
208
209#undef __FUNCT__
210#define __FUNCT__ "psi"
211
212/*++++++++++++++++++++++++++++++++++++++
213  This abstracts out the function for free energy
214  +latex+$\Psi(C_2,C_3)$,
215  +html+ Psi(<i>C_2</i>,<i>C_3</i>),
216  which is only used for integration/diagnostic purposes.
217
218  PetscScalar psi It returns the homogeneous free
219  energy.
220
221  PetscScalar C The
222  +latex+$C_2,C_3$
223  +html+ <i>C_2</i>, <i>C_3</i>
224  parameter it's a function of.
225
226  chtparm *thechternary chternary parameter structure.
227  ++++++++++++++++++++++++++++++++++++++*/
228
229static inline PetscScalar psi (PetscScalar C2,PetscScalar C3,chtparm *thechternary)
230{
231  int ierr;
232  /*    if (thechternary->polymer)
233        {
234        PetscScalar chi12=(0.218-0.276)/(1-0.622*C2/(1-C3));
235        PetscScalar chi13=3.5;
236        PetscScalar chi23=-1+0.5*C3;
237       
238        PetscScalar dchi12_dC3=(-0.058)/((1-0.622*C2/(1-C3))*(1-0.622*C2/(1-C3))) * 0.622*C2/((1-C3)*(1-C3));
239        PetscScalar dchi13_dC3=0.0;
240        PetscScalar dchi23_dC3=0.5;
241       
242        return -1.0/thechternary->m1-1.0/thechternary->m1*log(1.0-C2-C3)+1.0/thechternary->m3
243          +1.0/thechternary->m3*log(C3)-chi12*C2+chi23*C2+chi13*(1.0-2.0*C3-C2)+
244          (1-C2-C3)*C2*dchi12_dC3+C2*C3*dchi23_dC3+(1-C2-C3)*C3*dchi13_dC3;
245       
246          }*/
247       
248        if (thechternary->metal)
249          {
250            PetscScalar chiMM =3.3;
251            PetscScalar chiMC =2.5;
252            PetscScalar chipMM =0.0;
253            PetscScalar chipMC =0.0;
254
255            return C2*log(C2)+(1-C2)*log(1.0-C2)+C3*log(C3)+(1-C3)*log(1.0-C3) +
256              chiMC*C3*(1.0-C3)+  chiMM*C2*C3*(1.0-C2) - (C2-0.5)*(C3-0.5) +
257              chipMC*C3*C3*(1-C3)*(1-C3) + chipMM*C3*(C2*C2)*(1-C2)*(1-C2);
258          }
259        else
260        {
261                ierr=PetscPrintf(PETSC_COMM_WORLD,"Can't find the specified free energy formula!!!\n Please choose between 'polymer' and 'metal'\n");CHKERRQ (ierr);
262        }
263
264  return 0;
265}
266
267
268#undef __FUNCT__
269#define __FUNCT__ "psiprime22"
270
271/*++++++++++++++++++++++++++++++++++++++
272  This abstracts out the function for the second derivative of homogeneous free
273  energy with respect to
274  +latex+$C_2$: $\partial^2 G/\partial C_2^2 (C_2,C_3)$
275  +html+ <i>C</i><sub>2</sub>:
276  +html+ <i>d</i><sup>2</sup></i><i>G</i>/<i>dC</i><sub>2</sub><sup>2</sup>
277  +html+ (<i>C</i><sub>2</sub>, <i>C</i><sub>3</sub>)
278  so it can be easily modified.
279 
280  PetscScalar psiprime22 It returns the second derivative of homogeneous free
281  energy.
282
283  PetscScalar C2 Concentration of species 2.
284
285  PetscScalar C3 Concentration of species 3.
286
287  chtparm *thechternary chternary parameter structure.
288  ++++++++++++++++++++++++++++++++++++++*/
289
290static inline PetscScalar psiprime22 (PetscScalar C2, PetscScalar C3,
291                                     chtparm *thechternary)
292{
293  int ierr;
294
295  if (thechternary->polymer)
296    {
297      /*+ This is not ready yet. +*/
298
299      /* Polyurethane parameters */
300      PetscScalar chi_ns = 0.218 + 0.276/(1-0.622*C2/(1-C3));
301      PetscScalar chi_sp = 0.315;
302      PetscScalar chi_np = 3.53;
303
304      PetscScalar dchins_dC2 = 0.276/(1-0.622*C2/(1-C3))/(1-0.622*C2/(1-C3)) *
305        0.622*C2/(1-C3);
306      PetscScalar dchisp_dC2 = 0.0;
307      PetscScalar dchinp_dC2 = 0.0;
308
309      /*DPRINTF ("C2=%g, C3=%g, chi12=%g, chi13=%g, chi23=%g\n",
310               C2, C3, chi12,chi13,chi23);*/
311
312      return ((-1.0-pseudolog(1.-C2-C3))/thechternary->V_n +
313              (1.+pseudolog(C2))/thechternary->V_s +
314              (chi_ns*(1.-2.*C2-C3) +
315               dchins_dC2*(1.-C2-C3)*C2)/thechternary->V_n +
316              (chi_sp*C3 + dchisp_dC2*C2*C3)/thechternary->V_s +
317              (-chi_np*C3 + dchinp_dC2*(1.-C2-C3)*C3)/thechternary->V_n) *
318        thechternary->RT;
319    }
320
321  /* Metal electrolysis parameters and homogeneous free energy */
322  PetscScalar chiMM =3.3;
323  PetscScalar chiMC =2.5;
324  PetscScalar chipMM =0.0;
325  PetscScalar chipMC =0.0;
326
327  return log(C2) - log(1.0-C2) + chiMM*C3*(1.0-2.*C2)-(C3-0.5) + chipMM*C3*(2.*C2-6.*C2*C2 +4.*C2*C2*C2);
328}
329
330
331#undef __FUNCT__
332#define __FUNCT__ "chternary_first_setup"
333
334/*++++++++++++++++++++++++++++++++++++++
335  The basic setup, assigning the number of solved and temporary field
336  variables, the stencil width, and using options to set the parameters in the
337  +latex+{\tt chtparm}
338  +html+ <tt>chtparm</tt>
339  struct typedef.
340
341  int chternary_first_setup Returns zero (or an error code).
342
343  PetscTruth threedee Request support for 3-D.
344
345  int *vars Pointer to the number of solved field variables.
346
347  int *tempvars Pointer to the number of temporary field variables.
348
349  int *stencilwid Pointer to the stencil width.
350
351  AppCtx *data Pointer to the
352  +latex+{\tt AppCtx}
353  +html+ <tt>AppCtx</tt>
354  struct typedef, into whose
355  +latex+{\tt chtparm}
356  +html+ <tt>chtparm</tt>
357  structure this inserts parameters from the command line.
358  ++++++++++++++++++++++++++++++++++++++*/
359
360int chternary_first_setup
361(PetscTruth threedee, int *vars, int *tempvars, int *intvars, int *stencilwid, AppCtx *data)
362{
363  chtparm *thechternary = &data->thechternary;
364  int ierr;
365 
366  data->jacobian = PETSC_FALSE; /* Doesn't have one, but may at some point. */
367
368  thechternary->C2var  = (*vars)++;
369  thechternary->mu2var = (*tempvars)++;
370  thechternary->C3var  = (*vars)++;
371  thechternary->mu3var = (*tempvars)++;
372  thechternary->psivar = (*intvars)++;
373
374  ierr = PetscOptionsHasName(PETSC_NULL,"-cht_polymer",&thechternary->polymer);
375  CHKERRQ (ierr);
376  ierr = PetscOptionsHasName (PETSC_NULL, "-cht_metal", &thechternary->metal);
377  CHKERRQ (ierr);
378  if (!thechternary->polymer && !thechternary->metal)
379    {
380      ierr = PetscPrintf (PETSC_COMM_WORLD, "No chternary energy model selected, you need either -cht_polymer or -cht_metal\n");
381      return data->chternary = PETSC_FALSE;
382    }
383
384  /* Molar volumes in m^3, given for water-DMF-PVDF */
385  thechternary->V_n = 18e-6;
386  thechternary->V_s = 77.44e-6;
387  thechternary->V_p = 180e-6;
388  thechternary->RT  = 8.314*290.;
389  thechternary->ForceTerm=1.;
390  thechternary->uconst=0.;
391  thechternary->vconst=0.;
392  ierr = PetscOptionsGetScalar (PETSC_NULL, "-cht_Vn",
393                                &thechternary->V_n, PETSC_NULL); CHKERRQ (ierr);
394  ierr = PetscOptionsGetScalar (PETSC_NULL, "-cht_Vs",
395                                &thechternary->V_s, PETSC_NULL); CHKERRQ (ierr);
396  ierr = PetscOptionsGetScalar (PETSC_NULL, "-cht_Vp",
397                                &thechternary->V_p, PETSC_NULL); CHKERRQ (ierr);
398  ierr = PetscOptionsGetScalar (PETSC_NULL, "-cht_ForceTerm",
399                                &thechternary->ForceTerm, PETSC_NULL);
400  CHKERRQ (ierr);
401  ierr = PetscOptionsGetScalar (PETSC_NULL, "-cht_uUniform",
402                                &thechternary->uconst, PETSC_NULL);
403  CHKERRQ (ierr);
404  ierr = PetscOptionsGetScalar (PETSC_NULL, "-cht_vUniform",
405                                &thechternary->vconst, PETSC_NULL);
406  CHKERRQ (ierr);
407
408  /*+ This sets mobility and +latex+$K_{ij}$
409    +latex+{\tt chternary\_labels\_initcond()}.
410    +html+ <tt>chternary_labels_initcond()</tt>.
411    +*/
412
413  if (thechternary->polymer)
414    {
415      thechternary->mobility22 = 1e-8;
416      thechternary->mobility23 = 0;
417      thechternary->mobility32 = 0;
418      thechternary->mobility33 = 1e-9;
419      thechternary->K22= 1e-2;
420      thechternary->K23= 0;
421      thechternary->K32= 0;
422      thechternary->K33= 1e-2;
423    }
424  else
425    {
426      thechternary->mobility22 = 20;
427      thechternary->mobility23 = 0;
428      thechternary->mobility32 = 0;
429      thechternary->mobility33 = 200;
430      thechternary->K22= 1e-2;
431      thechternary->K23= 1e-2;
432      thechternary->K32= 1e-2;
433      thechternary->K33= 1e-2;
434    }
435
436  ierr = PetscOptionsGetScalar (PETSC_NULL,"-mobility_22",&thechternary->mobility22,PETSC_NULL); CHKERRQ (ierr);
437  ierr = PetscOptionsGetScalar (PETSC_NULL,"-mobility_23",&thechternary->mobility23,PETSC_NULL); CHKERRQ (ierr);
438  ierr = PetscOptionsGetScalar (PETSC_NULL,"-mobility_32",&thechternary->mobility32,PETSC_NULL); CHKERRQ (ierr);
439  ierr = PetscOptionsGetScalar (PETSC_NULL,"-mobility_33",&thechternary->mobility33,PETSC_NULL); CHKERRQ (ierr);
440
441  ierr = PetscOptionsGetScalar (PETSC_NULL, "-K_22", &thechternary->K22,PETSC_NULL); CHKERRQ (ierr);
442  ierr = PetscOptionsGetScalar (PETSC_NULL, "-K_23", &thechternary->K23,PETSC_NULL); CHKERRQ (ierr);
443  ierr = PetscOptionsGetScalar (PETSC_NULL, "-K_32", &thechternary->K32,PETSC_NULL); CHKERRQ (ierr);
444  ierr = PetscOptionsGetScalar (PETSC_NULL, "-K_33", &thechternary->K33,PETSC_NULL); CHKERRQ (ierr);
445  *stencilwid = PetscMax (*stencilwid, 2);
446
447  DPRINTF ("Vn=%g, Vs=%g, Vp=%g,\n", thechternary->V_n, thechternary->V_s,
448           thechternary->V_p);
449  DPRINTF ("mobility_22=%g, mobility_23=%g, mobility_32=%g, mobility_33=%g,\n",
450           thechternary->mobility22, thechternary->mobility23,
451           thechternary->mobility32, thechternary->mobility33);
452  DPRINTF ("K_22=%g, K_23=%g, K_32=%g, K_33=%g\n", thechternary->K22,
453           thechternary->K23, thechternary->K32, thechternary->K33);
454
455  /* Boundary condition parameters */
456  thechternary->hd_22=-1.;
457  thechternary->hd_33=-1.;
458  thechternary->c2ref=0.;
459  thechternary->c3ref=0.;
460  ierr = PetscOptionsGetScalar (PETSC_NULL,"-hd_22", &thechternary->hd_22,
461                                PETSC_NULL); CHKERRQ (ierr);
462  ierr = PetscOptionsGetScalar (PETSC_NULL,"-hd_33", &thechternary->hd_33,
463                                PETSC_NULL); CHKERRQ (ierr);
464  ierr = PetscOptionsGetScalar (PETSC_NULL,"-C2_ref", &thechternary->c2ref,
465                                PETSC_NULL); CHKERRQ (ierr);
466  ierr = PetscOptionsGetScalar (PETSC_NULL,"-C3_ref", &thechternary->c3ref,
467                                PETSC_NULL); CHKERRQ (ierr);
468
469  if (thechternary->hd_22 >= 0. || thechternary->hd_33 >= 0.)
470    DPRINTF ("BC parameters: hd_22=%g, hd_33=%g, C2_ref=%g, C3_ref=%g\n",
471             thechternary->hd_22, thechternary->hd_33, thechternary->c2ref,
472             thechternary->c3ref);
473
474  return 0;
475}
476
477
478#undef __FUNCT__
479#define __FUNCT__ "chternary_labels_initcond"
480
481/*++++++++++++++++++++++++++++++++++++++
482  This sets up the field variable labels, maximum stable explicit deltat, and
483  initial condition for the Cahn-Hilliard variables.
484
485  int chternary_labels_initcond Returns zero (or an error code).
486
487  PetscScalar *globalarray The global field array.
488
489  int nx Overall
490  +latex+$x$-width
491  +html+ <i>x</i>-width
492  of the global array.
493
494  int ny Overall
495  +latex+$y$-width
496  +html+ <i>y</i>-width
497  of the global array.
498
499  int nz Overall
500  +latex+$z$-width
501  +html+ <i>z</i>-width
502  of the global array.
503
504  int xm The
505  +latex+$x$-width
506  +html+ <i>x</i>-width
507  of the local part of the array.
508
509  int ym The
510  +latex+$y$-width
511  +html+ <i>y</i>-width
512  of the local part of the array.
513
514  int zm The
515  +latex+$z$-width
516  +html+ <i>z</i>-width
517  of the local part of the array.
518
519  int xs The (integer)
520  +latex+$x$-coordinate
521  +html+ <i>x</i>-coordinate
522  of the start of the local part of the array.
523
524  int ys The (integer)
525  +latex+$y$-coordinate
526  +html+ <i>y</i>-coordinate
527  of the start of the local part of the array.
528
529  int zs The (integer)
530  +latex+$z$-coordinate
531  +html+ <i>z</i>-coordinate
532  of the start of the local part of the array.
533
534  int vars Total number of field variables to be solved.
535
536  AppCtx *data Pointer to the
537  +latex+{\tt AppCtx}
538  +html+ <tt>AppCtx</tt>
539  struct typedef, whose
540  +latex+{\tt chtparm}
541  +html+ <tt>chtparm</tt>
542  structure this uses for various purposes.
543
544  PetscScalar *max_explicit_deltat Pointer to the largest allowable explicit
545  timestep size for this equation, which this function can set/modify.
546  ++++++++++++++++++++++++++++++++++++++*/
547
548int chternary_labels_initcond
549(PetscScalar *globalarray, int nx,int ny,int nz, int xm,int ym,int zm, int xs,
550 int ys,int zs, int vars, AppCtx *data, PetscScalar *max_explicit_deltat)
551{
552  chtparm *thechternary = &data->thechternary;
553  PetscScalar temp;
554  int ix,iy,iz, ierr;
555 
556  PetscScalar delta = PetscMin (data->xwid/nx, data->ywid/ny);
557    /* Set DA labels, constraints and symmetry types */
558  data->label[thechternary->C2var] = "C2";
559  data->label[thechternary->C3var] = "C3";
560
561  /*data->timestyle[thechternary->C2var] =  TIMESTEP_ONLY_VAR+thechternary->C2var;
562    data->timestyle[thechternary->C3var] = TIMESTEP_ONLY_VAR+thechternary->C3var; */
563  data->timestyle[thechternary->C2var]=TIME_CONST_BLEND_VAR+thechternary->C2var;
564  data->timestyle[thechternary->C3var] =TIME_CONST_BLEND_VAR+thechternary->C3var;
565
566  data->plot_types[thechternary->C2var] = FIELD_TERNARY;
567  data->plot_types[thechternary->C3var] = FIELD_TERNARY;
568  data->symmtypes[thechternary->C2var] =
569    SYMMETRY_MIRROR_PLANE * SYMMETRY_XMIN_START |
570    SYMMETRY_MIRROR_PLANE * SYMMETRY_YMIN_START |
571    SYMMETRY_MIRROR_PLANE * SYMMETRY_ZMIN_START |
572    SYMMETRY_MIRROR_PLANE * SYMMETRY_XMAX_START |
573    SYMMETRY_MIRROR_PLANE * SYMMETRY_YMAX_START |
574    SYMMETRY_MIRROR_PLANE * SYMMETRY_ZMAX_START;
575  data->symmtypes[thechternary->C3var] =
576    SYMMETRY_MIRROR_PLANE * SYMMETRY_XMIN_START |
577    SYMMETRY_MIRROR_PLANE * SYMMETRY_YMIN_START |
578    SYMMETRY_MIRROR_PLANE * SYMMETRY_ZMIN_START |
579    SYMMETRY_MIRROR_PLANE * SYMMETRY_XMAX_START |
580    SYMMETRY_MIRROR_PLANE * SYMMETRY_YMAX_START |
581    SYMMETRY_MIRROR_PLANE * SYMMETRY_ZMAX_START;
582
583  /* calculate explicit_deltat */
584  *max_explicit_deltat = PetscMin
585    (*max_explicit_deltat, .025/nx/nx/nx/thechternary->mobility22);
586  DPRINTF ("calculated deltat = %g\n",*max_explicit_deltat);
587
588  if (!data->load_data) /* Make sure we're not stepping on loaded data */
589    {
590      PetscTruth uniform, layers,trilayer, box,particles, triangle,sinusoidal;
591      int xmin,xmax, ymin,ymax, zmin,zmax, rank,size, random_counter;
592      PetscScalar fluct=0.005;
593      /* layer1=solution; layer2=bath */
594      PetscScalar C2_layer1=0.98,C3_layer1=0.98,C2_layer2=0.02,C3_layer2=0.02;
595      /* third layer for trilayer metal-electrolyte-metal structure */
596      PetscScalar C2_layer3=0.02, C3_layer3=0.98;
597
598      /*+ If we're not loading in data as the initial condition, the
599        +latex+$C2$ and $C3$
600        +html+ <i>C2</i> and <i>C3</i>
601        fields are initiated here.  There are several types of initial
602        conditions here which are chosen by command-line parameters:
603        +latex+\begin{itemize} \item
604        +html+ <ul><li>
605        +*/
606
607      ierr = PetscOptionsGetScalar (PETSC_NULL, "-cht_C2_layer1",
608                                    &C2_layer1, PETSC_NULL); CHKERRQ(ierr);
609      ierr = PetscOptionsGetScalar (PETSC_NULL, "-cht_C3_layer1",
610                                    &C3_layer1, PETSC_NULL); CHKERRQ(ierr);
611      ierr = PetscOptionsGetScalar (PETSC_NULL, "-cht_C2_layer2",
612                                    &C2_layer2, PETSC_NULL); CHKERRQ(ierr);
613      ierr = PetscOptionsGetScalar (PETSC_NULL, "-cht_C3_layer2",
614                                    &C3_layer2, PETSC_NULL); CHKERRQ(ierr);
615      ierr = PetscOptionsGetScalar (PETSC_NULL, "-cht_C2_layer3",
616                                    &C2_layer3, PETSC_NULL); CHKERRQ(ierr);
617      ierr = PetscOptionsGetScalar (PETSC_NULL, "-cht_C3_layer3",
618                                    &C3_layer3, PETSC_NULL); CHKERRQ(ierr);
619
620      /*+ Uniform initial condition is selected using the
621        +latex+{\tt -cht\_uniform\_C2} ard/or {\tt -cht\_uniform\_C3}
622        +html+ <tt>-cht_uniform_C2</tt> and/or <tt>-cht_uniform_C3</tt>
623        options.
624        +latex+\item
625        +html+ </li><li>
626        +*/
627      ierr = PetscOptionsHasName (PETSC_NULL,"-cht_uniform_C2",&uniform);
628      CHKERRQ (ierr);
629      if (uniform)
630        {
631          ierr = PetscOptionsGetScalar (PETSC_NULL, "-cht_uniform_C2",
632                                        &C2_layer2, PETSC_NULL); CHKERRQ(ierr);
633          DPRINTF ("Uniform C2=%g\n", C2_layer2);
634        }
635      ierr = PetscOptionsHasName (PETSC_NULL,"-cht_uniform_C3",&uniform);
636      CHKERRQ (ierr);
637      if (uniform)
638        {
639          ierr = PetscOptionsGetScalar (PETSC_NULL, "-cht_uniform_C3",
640                                        &C3_layer2, PETSC_NULL); CHKERRQ(ierr);
641          DPRINTF ("Uniform C3=%g\n", C3_layer2);
642        }
643
644      /* Uniform starting point: layer 2 composition */
645      for (iz=0; iz<(data->threedee ? zm: 1); iz++)
646        for (iy=0; iy<ym; iy++)
647          for (ix=0; ix<xm; ix++)
648            {
649              globalarray [((iz*ym + iy)*xm + ix)*vars + thechternary->C2var] =
650                C2_layer2;
651              globalarray [((iz*ym + iy)*xm + ix)*vars + thechternary->C3var] =
652                C3_layer2;
653            }
654
655      /*+ The
656        +latex+{\tt -cht\_particles}
657        +html+ <tt>-cht_particles</tt>
658        option is used for the ``colliding particles'' simulation.  This
659        creates a symmetric simulation with a square particle centered on the
660        middle of the
661        +latex+$x$-axis
662        +html+ <i>x</i>-axis
663        with width half that of the domain, which is to say, one-quarter of
664        that of the whole symmetric domain.
665        +latex+\item
666        +html+ </li><li>
667        +*/
668      ierr = PetscOptionsHasName (PETSC_NULL, "-cht_particles",&particles);
669      CHKERRQ (ierr);
670      if (particles)
671        {
672          xmin = nx/4;
673          xmax = 3*nx/4;
674          ymin = 0;
675          ymax = nx/4;
676          zmin = 0;
677          zmax = nz/4;
678          data->bcflags |= (XMIN_SYMMETRY | YMIN_SYMMETRY | ZMAX_SYMMETRY);
679        }
680
681      /*+ The "box" initial condition selected by
682        +latex+{\tt -cht\_box}
683        +html+ <tt>-cht_box</tt>
684        is a square half of the domain width, either in the center or, if the
685        +latex+$x$- and $y$-axes
686        +html+ <i>x</i>- and <i>y</i>-axes
687        are symmetry planes, then at the origin.
688        +latex+\item
689        +html+ </li><li>
690        +*/
691      ierr = PetscOptionsHasName (PETSC_NULL, "-cht_box", &box);
692      CHKERRQ (ierr);
693      if (box)
694        {
695          if ((data->bcflags & XMIN_SYMMETRY) &&
696              (data->bcflags & YMIN_SYMMETRY))
697            {
698              xmin = ymin = zmin = 0.;
699              xmax = nx/2;
700              ymax = ny/2;
701              zmax = nz/2;
702            }
703          else
704            {
705              xmin = nx/4;
706              xmax = 3*nx/4;
707              ymin = ny/4;
708              ymax = 3*ny/4;
709              zmin = nz/4;
710              zmax = 3*nz/4;
711            }
712        }
713
714      /*+ A two-layer initial condition in the
715        +latex+$y$-direction is chosen using the {\tt -cht\_layers}
716        +html+ <i>y</i>-direction is chosen using the <tt>-cht_layers</tt>
717        option, whose argument is the fraction of the domain which will be in
718        the ``bottom'' layer.
719        +latex+\item
720        +html+ </li><li>
721        +*/
722      ierr = PetscOptionsHasName (PETSC_NULL, "-cht_layers", &layers);
723      CHKERRQ (ierr);
724      if (layers)
725        {
726          PetscScalar layer_thickness;
727
728          ierr = PetscOptionsGetScalar (PETSC_NULL, "-cht_layers",
729                                        &layer_thickness, PETSC_NULL);
730          CHKERRQ (ierr);
731          DPRINTF ("Layer thickness: %g\n", layer_thickness);
732          xmin=ymin=zmin=0;
733          xmax=nx;
734          zmax=nz;
735          ymax=layer_thickness * ny;
736        }
737
738      /* Fill in the initial condition box */
739      for (iz = (data->threedee ? PetscMax (zs, zmin) : 0);
740           iz < (data->threedee ? PetscMin (zs+zm, zmax): 1); iz++)
741        for (iy=PetscMax(ys,ymin); iy<PetscMin(ys+ym,ymax); iy++)
742          for (ix=PetscMax(xs,xmin); ix<PetscMin(xs+xm,xmax); ix++)
743            {
744              globalarray [(((iz-zs)*ym + iy-ys)*xm + ix-xs)*vars +
745                           thechternary->C2var] =C2_layer1;
746              globalarray [(((iz-zs)*ym + iy-ys)*xm + ix-xs)*vars +
747                           thechternary->C3var] =C3_layer1;
748            }
749
750      /*+ A sinusoidal wave can crown the top of the bottom layer using option
751        +latex+{\tt -cht\_amplitude},
752        +html+ <tt>-cht_amplitude</tt>,
753        whose argument is the (absolute) amplitude of the sine wave with
754        wavelength equal to the domain width.
755        +latex+\item
756        +html+ </li><li>
757        +*/
758      ierr = PetscOptionsHasName (PETSC_NULL, "-cht_amplitude",&sinusoidal);
759      CHKERRQ (ierr);
760      if (sinusoidal && layers)
761        {
762          PetscScalar amplitude=delta, layer_thickness;
763          ierr = PetscOptionsGetScalar (PETSC_NULL, "-cht_amplitude",
764                                            &amplitude, PETSC_NULL);
765          ierr = PetscOptionsGetScalar (PETSC_NULL, "-cht_layers",
766                                        &layer_thickness, PETSC_NULL);
767
768          xmin=0;
769          xmax=nx;
770          zmin=0;
771          zmax=nz;           
772          ymin=0;
773          for (iz = (data->threedee ? PetscMax (zs, zmin) : 0);
774               iz < (data->threedee ? PetscMin (zs+zm, zmax): 1); iz++)
775            for (ix=PetscMax(xs,xmin); ix<PetscMin(xs+xm,xmax); ix++)
776              {
777                PetscScalar ymax_scalar = layer_thickness*ny +
778                  (amplitude/delta) *
779                  (2. + sin (2.*PETSC_PI*ix/nx) + sin (2.*PETSC_PI*iz/nz));
780                for (iy=PetscMax(ys,ymin); iy<PetscMin(ys+ym,ymax_scalar);
781                     iy++)
782                  {
783                    globalarray [(((iz-zs)*ym + iy-ys)*xm + ix-xs)*vars +
784                                 thechternary->C2var] = C2_layer1;
785                    globalarray [(((iz-zs)*ym + iy-ys)*xm + ix-xs)*vars +
786                                 thechternary->C3var] = C3_layer1;
787                    if (iy>ymax_scalar && iy !=ys)
788                      {
789                        globalarray [(((iz-zs)*ym + iy-1-ys)*xm + ix-xs)*vars +
790                                     thechternary->C2var] = ymax_scalar - (iy-1);
791                        globalarray [(((iz-zs)*ym + iy-1-ys)*xm + ix-xs)*vars +
792                                     thechternary->C3var] = ymax_scalar - (iy-1);
793                      }
794                  }
795              }
796        }
797
798      /*+ One can add a third layer using the
799        +latex+{\tt -cht\_trilayer}
800        +html+ <tt>-cht_trilayer</tt>
801        option, whose argument is the fraction of the domain above which will
802        reside the "top" layer.
803        +latex+\item
804        +html+ </li><li>
805        +*/
806      ierr = PetscOptionsHasName (PETSC_NULL, "-cht_trilayer", &trilayer);
807      CHKERRQ (ierr);
808      if (trilayer)
809        {
810          PetscScalar layer_thickness;
811
812          ierr = PetscOptionsGetScalar (PETSC_NULL, "-cht_trilayer",
813                                        &layer_thickness, PETSC_NULL);
814          CHKERRQ (ierr);
815          DPRINTF ("Layer thickness: %g\n", layer_thickness);
816
817          xmin = 0;
818          xmax = nx;
819          ymin = layer_thickness*ny;
820          ymax = ny;
821          zmin = 0;
822          zmax = nz;
823
824          for (iz = (data->threedee ? PetscMax (zs, zmin) : 0);
825               iz < (data->threedee ? PetscMin (zs+zm, zmax): 1); iz++)
826            for (iy=PetscMax(ys,ymin); iy<PetscMin(ys+ym,ymax); iy++)
827              for (ix=PetscMax(xs,xmin); ix<PetscMin(xs+xm,xmax); ix++)
828                {
829                  globalarray [(((iz-zs)*ym + iy-ys)*xm + ix-xs)*vars +
830                               thechternary->C2var] =C2_layer3;
831                  globalarray [(((iz-zs)*ym + iy-ys)*xm + ix-xs)*vars +
832                               thechternary->C3var] =C3_layer3;
833                }
834        }
835
836      /*+ If the option
837        +latex+{\tt -cht\_random\_fluct}
838        +html+ <tt>-cht_random_fluct</tt>
839        is selected without
840        +latex+{\tt -cht\_random\_center},
841        +html+ <tt>-cht_random_center</tt>,
842        then random fluctuations with uniform distribution of width twice the
843        fluctuation value are added to any other initial condition present.
844        +latex+\end{itemize}
845        +html+ </li></ul>
846        +*/
847
848      ierr = PetscOptionsGetScalar (PETSC_NULL, "-cht_random_fluct",&fluct,
849                                    PETSC_NULL); CHKERRQ (ierr);
850
851      /*+ Sometimes we only want to apply this fluctuation over the bottom
852        region of the simulation (e.g. for membranes), the option
853        +latex+{\tt -cht\_random\_layer}
854        +html+ <tt>-cht_random_layer</tt>
855        takes an argument from 0-1 and makes it operate over that fraction of
856        the domain. +*/
857
858      PetscScalar random_thickness=1.;
859      ierr = PetscOptionsGetScalar (PETSC_NULL, "-cht_random_layer",
860                                    &random_thickness, PETSC_NULL);
861      CHKERRQ (ierr);
862      ymax = ny * random_thickness;
863
864      /*+ For random fluctuations in this module, it's necessary to generate
865        different random sequences on each process, even though rand() should
866        return the same sequences.  So we create a random counter which takes
867        on values between zero and SMALL_PRIME-1; this is initialized to rank
868        times LARGE_PRIME (which should make it sort of random), and
869        incremented by MEDIUM_PRIME then re-moduloed to SMALL_PRIME for each
870        random number generated.  This counter, in turn, is divided by
871        SMALL_PRIME and the result added to rand()/RAND_MAX then modulo 1, in
872        order to give a "different random number" (at least at the resolution
873        of SMALL_PRIME) in each process.  This is not a great algorithm, but
874        seems to work, and in any case should do for the purpose needed here.
875
876        Eventually it might be good to make this resource available to the
877        whole code. +*/
878
879#define SMALL_PRIME 1571
880#define MEDIUM_PRIME 524287
881#define LARGE_PRIME 2147483647
882#define MY_RANDOM(pseudo_random) \
883  (((long long) (random_counter = (random_counter + MEDIUM_PRIME) % SMALL_PRIME) * RAND_MAX / SMALL_PRIME + pseudo_random) % RAND_MAX)
884      ierr = MPI_Comm_rank (PETSC_COMM_WORLD, &rank); CHKERRQ (ierr);
885      ierr = MPI_Comm_size (PETSC_COMM_WORLD, &size); CHKERRQ (ierr);
886      random_counter = ((long long) rank * LARGE_PRIME) % SMALL_PRIME;
887#ifdef DEBUG
888      ierr = PetscSynchronizedPrintf (PETSC_COMM_WORLD, "[%d] random_counter initialized to %d\n", rank, random_counter); CHKERRQ (ierr);
889      ierr = PetscSynchronizedFlush (PETSC_COMM_WORLD); CHKERRQ (ierr);
890#endif
891
892      ierr = PetscOptionsGetScalar (PETSC_NULL, "-cht_random_fluct",
893                                    &fluct, PETSC_NULL); CHKERRQ (ierr);
894      if (fluct != 0.)
895        for (iz = (data->threedee ? zs : 0);
896             iz < (data->threedee ? zs+zm: 1); iz++)
897          for (iy=ys; iy<PetscMin(ys+ym,ymax); iy++)
898            for (ix=xs; ix<xs+xm; ix++){
899              globalarray [(((iz-zs)*ym + iy-ys)*xm + ix-xs)*vars +
900                           thechternary->C2var] +=
901                -fluct + 2.*fluct*(PetscScalar)MY_RANDOM(rand())/RAND_MAX;
902              globalarray [(((iz-zs)*ym + iy-ys)*xm + ix-xs)*vars +
903                           thechternary->C3var] +=
904                -fluct + 2.*fluct*(PetscScalar)MY_RANDOM(rand())/RAND_MAX;
905            }
906    }
907
908  return 0;
909}
910
911#define C2(point) (x [vars*(point) + thechternary->C2var])
912#define C3(point) (x [vars*(point) + thechternary->C3var])
913#define C2func(point) (func [vars*(point) + thechternary->C2var])
914#define C3func(point) (func [vars*(point) + thechternary->C3var])
915#define mu2(point) (temp [tempvars*(point) + thechternary->mu2var])
916#define mu3(point) (temp [tempvars*(point) + thechternary->mu3var])
917#define Psi(point) (integ [intvars*(point) + thechternary->psivar])
918
919#define u(point) (x [vars*(point) + thevortex->uvar])
920#define v(point) (x [vars*(point) + thevortex->vvar])
921
922#define sigma(point) (temp[tempvars*(point) + thepotential->sigmavar])
923#define V(point) (x[vars*(point) + thepotential->Vvar])
924
925#define pu(point) (x [vars*(point) + thepressure->uvar])
926#define pv(point) (x [vars*(point) + thepressure->vvar])
927#define pw(point) (x [vars*(point) + thepressure->wvar])
928
929#undef __FUNCT__
930#define __FUNCT__ "chternary_integrate_interior_line"
931
932/*++++++++++++++++++++++++++++++++++++++
933  This calculates the chternary integration variable
934  +latex+$\Psi$,
935  +html+ <i>Psi</i>,
936  which is the free energy.
937 
938  int chternary_integrate_interior_line Returns zero (or an error code).
939
940  PetscScalar *x Array with the "real" field variables.
941
942  PetscScalar *integ Array with the integration variables.
943
944  int points Number of points at which to calculate the integration variables.