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

KineticLinearSolverCustom.groovy

import edu.gcsc.vrl.ug.api.*;
import eu.mihosoft.vrl.annotation.*;
import eu.mihosoft.vrl.math.Trajectory;
import eu.mihosoft.vrl.types.CanvasRequest;

import java.io.File;
import java.io.FileFilter;

/**
 * Script to perform a Kinetic Problem.
 *
 * Author: M. Hoffer, A. Vogel
 */
@ComponentInfo(name="KineticLinearSolver (custom version)", category="Custom")
public class KineticLinearSolverCustom 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="Data Evaluation|false")
            @ParamInfo(name="Volume Integration Subset  ", style="ugx-subset-selection-array", options="ugx_globalTag=\"TheFile\"; minArraySize=0")
                    String[] vEvalSubset,
            @ParamGroupInfo(group="Data Evaluation")
            @ParamInfo(name="Flux Integration Subset  ", style="ugx-subset-selection-array", options="ugx_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("Geometric MultiGrid".equals(precondName)){
            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     ("Implicit Euler".equals(timeSolverName)) timeDisc.set_theta(1.0);
        else if("Crank-Nicolson".equals(timeSolverName)) timeDisc.set_theta(0.5);
        else if("Explicit Euler".equals(timeSolverName)) 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);
        I_AssembledLinearOperator A = new AssembledLinearOperator(timeDisc);
        I_GridFunction b = 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     ("BiCGStab".equals(solverName)) solver = new BiCGStab();
        else if("Conjugate Gradient".equals(solverName)) solver = new CG();
        else if("Linear Solver".equals(solverName)) solver = new LinearSolver();
        else if("LU".equals(solverName)) solver = new LU();
        else errorExit("Cannot find linear solver: " + solverName);
        solver.set_convergence_check(convCheck);

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

            if     ("Jacobi".equals(precondName)) {
                precond = new Jacobi();
            }
            else if("Gauss-Seidel".equals(precondName)) {
                precond = new GaussSeidel();
            }
            else if("Incomplete LU".equals(precondName)) {
                precond = new ILU();
            }
            else if("Geometric MultiGrid".equals(precondName)) {
                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);
            }

            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 !?");
            }

        }

        //--------------------------------------------------------------------------------
        //--  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");
        }
        for(int fct = 0; fct <  approxSpace.const__num_fct(); ++fct){
            F_Interpolate.invoke(StartValues[fct], u, approxSpace.const__name(fct), time);
        }

        //-- 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);

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

        // 1. assemble matrix and rhs
        timeDisc.assemble_linear(A, b);

        // 3. init solver for linear Operator
        solver.init(A, u);

        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 + 1e-12 < 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 + 1e-12) {
                    dt = timeEnd - time + 1e-20;
                }

                // 1. assemble matrix and rhs
                timeDisc.prepare_step(solTimeSeries, dt);
                timeDisc.assemble_rhs(b);

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

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

                    timeDisc.prepare_step(solTimeSeries, dt);
                    timeDisc.assemble_linear(A, b);
                    solver.init(A, u);

                }else{
                    bSuccess = true;
                }
            }

            if(dt < maxStepSize){
                timeDisc.prepare_step(solTimeSeries, maxStepSize);
                timeDisc.assemble_linear(A, b);
                solver.init(A, u);
            }

            //-- 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 KineticLinearSolver: ", 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

}