KineticSolver-UG4 component for VRL by A.Vogel and M.Hoffer

KineticSolver.groovy

/**
 * Script to perform a Kinetic Problem
 * 
 * Author: M. Hoffer, A. Vogel
 */
import eu.mihosoft.vrl.annotation.*;
import eu.mihosoft.vrl.types.*;
import edu.gcsc.vrl.ug.api.I_ApproximationSpace;
import edu.gcsc.vrl.ug.api.I_UserNumber;
import edu.gcsc.vrl.ug.api.F_Integral;
import edu.gcsc.vrl.ug.api.*;
import eu.mihosoft.vrl.math.Trajectory;
import java.io.File;
import java.io.FileFilter;

@ComponentInfo(name = "KineticSolver (custom version)", category = "Custom")
public class KineticSolverCustom implements java.io.Serializable {

    private static final long serialVersionUID = 1L;
    private boolean stopSolver;
    private Trajectory[] vEvalTrajectory;
    private transient CanvasRequest cReq;

    @MethodInfo(valueStyle = "multi-out", interactive = true)
    @OutputInfo(
            style = "multi-out",
            elemNames = ["VTU-File", "Trajectories"],
            elemTypes = [File.class, Trajectory[].class]
    )
    public Object[] invoke(
            CanvasRequest ThecReq,
            @ParamGroupInfo(group = "Problem Setup")
            @ParamInfo(name = "Problem-Setup", style = "default") I_DomainDiscretization domainDisc,
            @ParamGroupInfo(group = "Problem Setup")
            @ParamInfo(name = "Approximation-Space", style = "default") I_ApproximationSpace approxSpace,
            @ParamGroupInfo(group = "Domain")
            @ParamInfo(name = "Output Filename:", style = "save-dialog", options = "endings=[\"vtu\"]") java.io.File fileOut,
            @ParamGroupInfo(group = "Domain")
            @ParamInfo(name = "Erase old Files", options = "value=true") boolean bEraseOldFiles,
            @ParamGroupInfo(group = "Domain")
            @ParamInfo(name = "Output-Data", style = "default") I_VTKOutput out,
            @ParamGroupInfo(group = "Domain")
            @ParamInfo(name = "Number Refinements", style = "default", options = "value=3") int numRefs,
            @ParamGroupInfo(group = "Time|true|Time Solver Parameters")
            @ParamInfo(name = "Time Solver", style = "selection", options = "value=[\"Implicit Euler\", \"Crank-Nicolson\", \"Explicit Euler\"]") String timeSolverName,
            @ParamGroupInfo(group = "Time")
            @ParamInfo(name = "Start Time", style = "default", options = "value=0.00D") double timeStart,
            @ParamGroupInfo(group = "Time")
            @ParamInfo(name = "End Time", style = "default", options = "value=0.1D") double timeEnd,
            @ParamGroupInfo(group = "Time")
            @ParamInfo(name = "Maximal Time Step Size", style = "default", options = "value=0.01D") double maxStepSize,
            @ParamGroupInfo(group = "Time")
            @ParamInfo(name = "Minimal Time Step Size", style = "default", options = "value=0.0001D") double minStepSize,
            @ParamGroupInfo(group = "Time")
            @ParamInfo(name = "Time Step Reduction Factor", style = "default", options = "value=0.25D") double stepReductionFactor,
            @ParamGroupInfo(group = "Time")
            @ParamInfo(name = "Start Values", style = "default") I_UserNumber[] StartValues,
            @ParamGroupInfo(group = "Linear Solver Setup|true|Solver Parameters")
            @ParamInfo(name = "Linear Solver", style = "selection", options = "value=[\"Linear Solver\", \"Conjugate Gradient\", \"BiCGStab\", \"LU\"]") String solverName,
            @ParamGroupInfo(group = "Linear Solver Setup")
            @ParamInfo(name = "Preconditioner", style = "selection", options = "value=[\"Geometric MultiGrid\", \"Jacobi\", \"Gauss-Seidel\", \"Incomplete LU\", \"None\"]") String precondName,
            @ParamGroupInfo(group = "Linear Solver Setup")
            @ParamInfo(name = "Maximal Number Iterations", style = "default", options = "value=100") int maxNumIter,
            @ParamGroupInfo(group = "Linear Solver Setup")
            @ParamInfo(name = "Minimal Residuum Norm", style = "default", options = "value=1E-10") double absTol,
            @ParamGroupInfo(group = "Nonlinear Solver Setup")
            @ParamInfo(name = "Maximal Number Iterations", style = "default", options = "value=20") int maxNumIterNonLinear,
            @ParamGroupInfo(group = "Nonlinear Solver Setup")
            @ParamInfo(name = "Minimal Residuum Norm", style = "default", options = "value=1E-8D") double absTolNonLinear,
            @ParamGroupInfo(group = "Data Evaluation|false")
            @ParamInfo(name = "Integration Subset  ", style = "ugx-subset-selection-array", options = "globalTag=\"TheFile\"; minArraySize=0") String[] vEvalSubset,
            @ParamGroupInfo(group = "Data Evaluation")
            @ParamInfo(name = "Flux Integration Subset  ", style = "ugx-subset-selection-array", options = "globalTag=\"TheFile\"; minArraySize=0") String[] vEvalFluxSubset) {
        if (bEraseOldFiles) {
            eraseAllFilesInFolder(fileOut, "vtu");
        }

        cReq = ThecReq;

        String fileNameOut = fileOut.getAbsoluteFile().getAbsolutePath();

        // set abortion flag to false initially (can be changed using stopSolver-Method)
        stopSolver = false;

        I_Domain dom = approxSpace.domain();
        I_IRefiner refiner = F_GlobalDomainRefiner.invoke(dom);
        for (int i = 0; i < numRefs; ++i) {
            refiner.refine();
        }

        if (precondName.equals("Geometric MultiGrid")) {
            approxSpace.init_levels();
        }
        approxSpace.init_top_surface();
        approxSpace.const__print_statistic();

        F_OrderCuthillMcKee.invoke(approxSpace, true);

        //-- Create Time Discretization
        I_ThetaTimeStep timeDisc = new ThetaTimeStep(domainDisc);
        if (timeSolverName.equals("Implicit Euler")) {
            timeDisc.set_theta(1.0);
        } else if (timeSolverName.equals("Crank-Nicolson")) {
            timeDisc.set_theta(0.5);
        } else if (timeSolverName.equals("Explicit Euler")) {
            timeDisc.set_theta(0.0);
        } else {
            errorExit("Cannot find time solver: " + timeSolverName);
        }

 //        timeDisc.set_theta(1.0); //-- 1.0 is implicit euler
        //-- create operator from discretization
        I_AssembledOperator op = new AssembledOperator();
        op.set_discretization(timeDisc);
        op.init();

        I_MGSubsetHandler sh = dom.subset_handler();

        //--------------------------------------------------------------------------------
        //--  Algebra
        //--------------------------------------------------------------------------------
        //-- matrix and vectors
        I_GridFunction u = new GridFunction(approxSpace);

        //-- create Convergence Check
        I_ConvCheck convCheck = new ConvCheck();
        convCheck.set_maximum_steps(maxNumIter);
        convCheck.set_minimum_defect(absTol);
        convCheck.set_reduction(1e-90);

        //-- create Solver with ilu preconditioner
        I_ILinearOperatorInverse solver = null;
        if (solverName.equals("BiCGStab")) {
            solver = new BiCGStab();
        } else if (solverName.equals("Conjugate Gradient")) {
            solver = new CG();
        } else if (solverName .equals( "Linear Solver")) {
            solver = new LinearSolver();
        } else if (solverName .equals( "LU")) {
            solver = new LU();
        } else {
            errorExit("Cannot find solver: " + solverName);
        }
        solver.set_convergence_check(convCheck);

        // -- create GMG Preconditioner
        if (!precondName.equals("None") && !solverName.equals("LU")) {
            I_ILinearIterator precond = null;

            if (precondName .equals( "Jacobi")) {
                precond = new Jacobi();
            } else if (precondName .equals( "Gauss-Seidel")) {
                precond = new GaussSeidel();
            } else if (precondName .equals( "Incomplete LU")) {
                precond = new ILU();
            } else if (precondName .equals( "Geometric MultiGrid")) {
                precond = new GeometricMultiGrid(approxSpace);
                ((GeometricMultiGrid) precond).set_discretization(timeDisc);
                ((GeometricMultiGrid) precond).set_base_level(0);
                ((GeometricMultiGrid) precond).set_base_solver(new LU());
                ((GeometricMultiGrid) precond).set_smoother(new ILU());
                ((GeometricMultiGrid) precond).set_cycle_type(1);
                ((GeometricMultiGrid) precond).set_num_presmooth(3);
                ((GeometricMultiGrid) precond).set_num_postsmooth(3);
            } else {
                errorExit("Cannot find preconditioner: " + precondName);
            }

//            solver.set_preconditioner(precond);
            try {
                ((I_IPreconditionedLinearOperatorInverse) solver).set_preconditioner(precond);
            } catch (Exception e) {
                // in case of solver is LU
                // there will be (maybe) a CastExeption
                System.err.println(getClass().getSimpleName() + ": solver is LU therefore CastException !?");
            }
        }

        I_ConvCheck newtonConvCheck = new ConvCheck();
        newtonConvCheck.set_maximum_steps(maxNumIterNonLinear);
        newtonConvCheck.set_minimum_defect(absTolNonLinear);
        newtonConvCheck.set_reduction(1e-20);
        newtonConvCheck.set_verbose(true);

        //-- create Newton Solver
        I_NewtonSolver newtonSolver = new NewtonSolver();
        newtonSolver.set_linear_solver(solver);
        newtonSolver.set_convergence_check(newtonConvCheck);
        newtonSolver.init(op);

        //--------------------------------------------------------------------------------
        //--  Output
        //--------------------------------------------------------------------------------
        // this is different from the original script:
        // we remove .vtu because ug's vtk function does not allow
        // filename extensions
        if (fileNameOut.endsWith(".vtu")) {
            fileNameOut = fileNameOut.substring(
                    0, fileNameOut.length() - 4);
        }

        //-- start
        double time = timeStart;
        int step = 0;

        //-- set initial value
        F_Print.invoke("Interpolation start values.\n");
        if (approxSpace.const__num_fct() > StartValues.length) {
            errorExit("Wrong number of start values passed");
        }
        
        System.out.println("approxSpace.const__num_fct() = "+ approxSpace.const__num_fct());
        
        for (int fct = 0; fct < approxSpace.const__num_fct(); ++fct) {
            F_Interpolate.invoke(StartValues[fct], u, approxSpace.const__name(fct), time);
            
//            //DEBUG START
//            //writes a file in the "studio-bundle folder"/Resources/Java/"filename.vec"
//            F_SaveVectorForConnectionViewer.invoke(u, "u_interpolated.vec");
//            
//            System.out.println("StartValues[fct], u, approxSpace.const__name(fct), time");
////            ConstUserNumber2d  n = (ConstUserNumber2d) StartValues[fct];
////            n.const__print();
////            PrintUserNumber2d pun = new PrintUserNumber2d();
////            F_PrintUserDataValue.invoke(StartValues[fct], , time, 0);
//            System.out.println(StartValues[fct].toString() +","+ u+","+ approxSpace.const__name(fct)+","+ time);
//            
//            if(StartValues[fct] instanceof VRLUserNumber2d){
//                VRLUserNumber2d number = (VRLUserNumber2d) StartValues[fct];
////                F_Print
//            }else if(StartValues[fct] instanceof ConstUserNumber2d){
//                ConstUserNumber2d number = (ConstUserNumber2d) StartValues[fct];
//                System.out.println("const_print = ");
//                number.const__print();
//            }
            //            //DEBUG END       
        }
        
        //-- write start solution
        F_Print.invoke("Writing start values.\n");
        out.print(fileNameOut, u, step, time);

        //-- create new grid function for old value
//        I_GridFunction uOld = u.clone();
        I_GridFunction uOld = u.const__clone();

        //-- store grid function in vector of  old solutions
        I_SolutionTimeSeries solTimeSeries = new SolutionTimeSeries();
        solTimeSeries.push(uOld, time);

        timeDisc.prepare_step(solTimeSeries, maxStepSize);

        vEvalTrajectory = new Trajectory[vEvalSubset.length + vEvalFluxSubset.length];
        for (int i = 0; i < vEvalSubset.length; ++i) {
            vEvalTrajectory[i] = new Trajectory(""+vEvalSubset[i]);
            double val = F_Integral.invoke(u, approxSpace.const__name(0), vEvalSubset[i]);
            F_Print.invoke("Volume Integral on Subset " + vEvalSubset[i] + ":" + val + "\n");
            vEvalTrajectory[i].add(time, val);
        }

        for (int i = 0; i < vEvalFluxSubset.length; ++i) {
            vEvalTrajectory[i + vEvalSubset.length] = new Trajectory(""+vEvalFluxSubset[i]);
            double val = -1.0 * F_IntegrateDiscFlux.invoke(timeDisc, u, approxSpace.const__name(0), vEvalFluxSubset[i]);
            F_Print.invoke("Flux Integral on Subset " + vEvalFluxSubset[i] + ":" + val + "\n");
            vEvalTrajectory[i + vEvalSubset.length].add(time, val);
        }

        while (time < timeEnd) {
            step += 1;
            F_Print.invoke("++++++ TIMESTEP " + step + "  BEGIN, time: " + time + " ++++++\n");

            //-- choose time step
            double dt = maxStepSize;

            boolean bSuccess = false;

            while (bSuccess == false) {

                if (time + dt > timeEnd) {
                    dt = timeEnd - time + 1e-20;
                }

                //-- setup time Disc for old solutions and timestep
                timeDisc.prepare_step(solTimeSeries, dt);

                //-- prepare newton solver
                if (newtonSolver.prepare(u) == false) {
                    F_Print.invoke("Newton solver failed at step " + step + ".");
                    errorExit("Newton Solver: Solving failed.");
                }

                //-- apply newton solver
                if (newtonSolver.apply(u) == false) {
                    dt = dt * stepReductionFactor;
                    F_Print.invoke("\n++++++ Newton solver failed. Trying decreased stepsize " + dt);

                    if (dt < minStepSize) {
                        F_Print.invoke("++++++ Time Step to small. Cannot solve problem.");
                        errorExit("Newton Solver: Solving failed.");
                    }

                } else {
                    bSuccess = true;
                }
            }

            //-- update new time
            time = solTimeSeries.const__time(0) + dt;

            for (int i = 0; i < vEvalSubset.length; ++i) {
                double val = F_Integral.invoke(u, approxSpace.const__name(0), vEvalSubset[i]);
                F_Print.invoke("Integral on Subset " + vEvalSubset[i] + ":" + val + "\n");
                vEvalTrajectory[i].add(time, val);
            }
            for (int i = 0; i < vEvalFluxSubset.length; ++i) {
                double val = -1.0 * F_IntegrateDiscFlux.invoke(timeDisc, u, approxSpace.const__name(0), vEvalFluxSubset[i]);
                F_Print.invoke("Flux Integral on Subset " + vEvalFluxSubset[i] + ":" + val + "\n");
                vEvalTrajectory[i + vEvalSubset.length].add(time, val);
            }

            //-- plot solution
            out.print(fileNameOut, u, step, time);

            //-- get oldest solution
            I_Vector oldestSol = solTimeSeries.oldest();

            //-- copy values into oldest solution (we reuse the memory here)
            F_VecScaleAssign.invoke(oldestSol, 1.0, u);

            //-- push oldest solutions with new values to front, oldest sol pointer is poped from end
            solTimeSeries.push_discard_oldest(oldestSol, time);

            F_Print.invoke("++++++ TIMESTEP " + step + "  END, time: " + time + " ++++++\n");

            if (stopSolver) {
                F_Print.invoke("\n -------- SOLVER STOPPED --------\n");
                break;
            }
        }

        //-- end timeseries, produce gathering file
        out.write_time_pvd(fileNameOut, u);

        return [new File(fileNameOut + ".vtu"), vEvalTrajectory];
    }

    public void stopSolver() {
        stopSolver = true;
    }

    private void errorExit(String s) {
        eu.mihosoft.vrl.system.VMessage.exception("Setup Error in KineticSolver: ", s);
    }

    @MethodInfo(name = "", valueName = "Trajectories", valueStyle = "default", valueOptions = "", interactive = false)
    public Trajectory[] getTrajectories() {
        return vEvalTrajectory;
    }

    static private void eraseAllFilesInFolder(
            File file,
            final String ending) {

//      int dot = file.getPath().lastIndexOf(".");
//      int sep = file.getPath().lastIndexOf(File.separator);
//      final String fileNameStartsWith = file.getPath().substring(sep + 1, dot);
        File folder = file.getParentFile();

        String path = file.getAbsolutePath();
        String fileNameStartsWith = path;
        if(path.toLowerCase().endsWith(ending.toLowerCase())) {
            fileNameStartsWith = path.substring(0,path.length()-ending.length());
        }

        // check if filename given
        if(fileNameStartsWith.isEmpty()) return;

        // we need to declare a final copy of fileNameStartsWith to make it accessible to
        // the inner filter class
        final String fileNameStartsWithFinal = fileNameStartsWith;
        // check that folder exists
        if (folder != null && folder.isDirectory()) {
            for (File f : folder.listFiles(new FileFilter() {

                @Override
                public boolean accept(File pathName) {
                    boolean fileAccept = pathName.getName().toLowerCase().endsWith(ending);

                    //int dot1 = pathName.getPath().lastIndexOf(".");
                    //int sep1 = pathName.getPath().lastIndexOf(File.separator);
                    String fileName = pathName.getPath().substring(0,pathName.getPath().length()-ending.length());

                    boolean nameAccept = fileName.startsWith(fileNameStartsWithFinal + "_t");

                    return fileAccept && nameAccept;
                }
            }))
            { // start for
                if (f.isFile() && f.exists()) { // start if exists
                    try{
                        //System.out.println("DELETING FILE: " + f.getAbsoluteFile());
                        f.delete();
                    } catch (Exception ex) {
                        // ... do nothing ... ignore ...
                    }
                } // end if exists
            } // end for
        } // end if folder

    } // end method
}