dmorse / pscfpp

Polymer Self-Consistent Field Theory (C++/CUDA version)
https://pscf-home.cems.umn.edu
GNU General Public License v3.0
25 stars 20 forks source link

Having a hard time reusing output files from previous simulations #153

Closed sourproton closed 7 months ago

sourproton commented 7 months ago

Hello, I'm here seeking help, and would like to share some difficulty I've been encountering when trying to reuse output files from previous simulations.

In fact, I've been using pscfpp for a short time, coming from the Fortran version of the software and having recently switched to experiment with some of the new features.

One thing I used to do with the Fortran version is to generate a density field file with one of the examples, then tune the params file to match the system I wanted to simulate and pass the density field from the previous simulation as an initial guess for the new simulation. The quasi-Newton-Raphson iteration algorithm was robust enough to converge, even if severe parameter changes were made, for example changing the chi value from 15 to 50 in a diblock, f=0.5 diblock polymer melt.

I'm having a hard trying trying to replicate the same methodology with pscf_pc1. From the example examples/pc/diblock/lam, I can generate symmetry-adapted files of the density and the chemical potential fields, c.bf and w.bf respectively, and pass them as inputs, either directly with w.bf or with ESTIMATE_W_FROM_C and c.bf, and redo the simulation without changing the parameters. However, a simple change in chi from 12 to 15 yields a homogeneous solution, opposite from the expected which is a more segregated density field when compared to the previous solution with chi = 12.

Another thing I tried is to keep the initial guess from the example, and change fA=0.56 and fB=0.44 to fA=fB=0.5. The solution converges, but when I try to pass this converged solution as the inital guess for a second simulation which changes the value of chi, the algorithm diverges and I can't get a solution.

The documentation states "the initial guess for the volume fraction fields does not need to be perfect, just close enough to the final structure to allow a chosen iteration algorithm to converge". Ultimately, I would also like to enter "real" values in the parameter file, as opposed to using nondimensionalized units. I would like to know if I'm overestimating the robustness of the Anderson-mixing iteration algorithm, or if I'm having the wrong approach when trying to modify the parameters and still have a converged solution.

Thank you for the active development of PSCF,

Mateus

sourproton commented 7 months ago

Here is a step-by-step description of what I'm trying:

  1. Run the example as is cd examples/pc/diblock/lam ./run -> this runs the simulation with the initial guess and parameters provided by the example

  2. Reuse the last field as initial guess cp out/w.bf in/w1.bf modify commands to read in/w1.bf as the initial guess ./run -> this runs the same simulation, but with the last density field as an initial guess. Nothing changed and as expected the solution converges

  3. Reuse the last field as initial guess and change fA and fB cp out/w.bf in/w2.bf modify commands to read in/w2.bf as the initial guess modify params's blocks from

    0 0.56
    1 0.44

    to

    0 0.5
    1 0.5

    ./run -> this runs the simulation, but with the last density field as an initial guess and the new values of fA and fB. This converges too

  4. Reuse the last field as initial guess and change chi from 12 to any value bigger than 12.005 cp out/w.bf in/w3.bf modify commands to read in/w3.bf as the initial guess modify chi from 1 0 12.0 to 1 0 13.0 ./run -> this runs the simulation, but with the last density field as an initial guess and the new value of chi.

The previous attempt fails to converge in 100 iterations, which is the current value of maxItr. Changing it from 100 to 1000 makes the solution converge to an homogeneous density field where every point has phiA and phiB = 0.5.

If instead of modifying maxItr, we modify maxHist from 10 to 50, the error converges to nan on iteration 14.

dmorse commented 7 months ago

Suggestion: Have you tried using the SWEEP functionality to change parameters in smaller steps, rather than trying to make such a drastic change in parameters in one step? That is what I would do if I wanted to change parameters and faced problems with convergence when I tried to do it in one step by the procedure you describe. Please let me know if this helps, or if there is some reason you can't or don't want to use a sweep.

Comment: I'm not surprised that you found the old Fortran quasi-Newton-Raphson iteration algorithm more robust than any iterators used in the current version. A lot of work was put into that algorithm (mostly by Jian Qin) in order to make it robust, and I don't think we haven't payed the same sort of attention to robustness in designing the newer Anderson-mixing algorithms. We have not ported the quasi-Newton-Raphson algorithm from Fortran to C++ primarily because it uses too much memory to be applied to large 3D problems, and because we've found the Anderson mixing algorithms which use less memory to be good substitutes for more challenging applications involving 3D structures. The old Fortran NR algorithm may, however, have been nearly ideal for problems involving 1D lamellar phases of the type used in your example,. I'm open to considering or suggesting ideas for ways to modify existing algorithms to make them more robust, or to accepting pull request from any user who wants to take the time to port the quasi-Newton algorithm from Fortran to C++.

sourproton commented 7 months ago

Thank you for the quick response. Indeed I thought about the SWEEP functionality. The main reason I didn't use it yet is because I don't master its syntax yet, but I'm aware that slowly incrementing the variables I want to explore while inter/extrapolating the initial guesses based on the last results could be interesting.

I think one of the problems of my approach described above was to reuse the chemical potential field from one simulation as the initial guess of another simulation with different chis. Since then I've had some success by using the density field from the previous simulation and using the function ESTIMATE_W_FROM_C to have the new simulation estimate the new initial guess of the chemical potential field with the new values of chi. I've noticed the convergence is greatly influenced by the initial guess of the cell parameter, but with a bit of trial and error based on the radius of gyration of the chain, it's possible to find a cell parameter that converges to a segregated solution.

I'll close this issue for now but will keep in mind that you are open to ideas and pull requests. Thank you.